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Abstract 

Helicopter  flight,  relying  on  rotary  motion  of  a  complex  mechanical  system,  is 
predisposed  to  vibration.  While  the  helicopter  has  many  sources  of  vibrations 
the  main  rotor  system  generates  by  far  the  largest  magnitude  vibrations,  and  can 
render  the  vehicle  inoperable  if  left  unaddressed.  Thus,  proper  maintenance  to  reduce 
vibrations  is  essential  to  the  safe  operation  of  any  helicopter.  This  maintenance, 
however,  is  costly  and  time  consuming.  Improving  the  maintenance  procedure  for 
balancing  the  main  rotor  system  has  been  an  area  of  active  interest  since  the  inception 
of  the  helicopter.  However,  the  state  of  the  art  in  rotor  balancing  still  requires  several 
iterations  of  rotor  adjustments,  each  necessitating  a  separate  test  flight  and  then  time 
consuming  maintenance,  to  reduce  the  vibrational  level  to  an  acceptable  amount.  This 
research  provides  the  basis  for  an  improved  rotor  vibrational  reduction  methodology 
that  significantly  reduces  the  number  of  adjustment  iterations  required  to  reduce  main 
rotor  vibrations. 

To  address  these  issues,  it  was  the  intent  of  this  research  to  develop  an  on-line, 
linear  time  periodic  rotor  vibration  controller.  The  Cramer-Rao  bound  was  developed 
for  a  linear  time  periodic  system  in  order  to  identify  the  quality  of  identified  system 
parameters  that  are  used  in  system  models  for  controller  development.  The  methods 
developed  in  this  work  have  allowed  model  parameters  to  be  verified  for  accuracy  and 
likewise  adjusted  to  improve  controller  accuracy.  To  achieve  these  goals  several  steps 
were  undertaken  as  enumerated  below. 

1.  Describe  the  methodology  defined  by  Wereley  [42]  to  model  a  linear  system  with 
time  periodic  coefficients  as  a  state  space  model,  in  a  manner  similar  to  a  linear 
time  invariant  system. 


IV 


2.  Develop  the  Cramer- Rao  Bound  to  validate  the  parameters  of  the  linear  time 
parametric  system  in  state  space  form,  as  in  the  case  of  a  helicopter  rotor  model. 

3.  Model  a  helicopter  rotor  system  which  incorporates  time  periodic  system  coef¬ 
ficients  to  accurately  describe  the  system  in  forward  flight. 

4.  Using  the  linear  time  periodic  state  space  model,  perform  system  identification 
of  the  main  rotor  system  to  identify  the  time  periodic  parameters  of  the  model. 

5.  Develop  an  optimal  control  methodology  for  a  linear  time  periodic  rotor  model 
as  to  provide  a  vibration  smoothing  solution  for  an  unbalanced  model. 
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Rotorcraft  Smoothing  via  Linear  Time  Periodic  Methods 


I.  Introduction 

Helicopter  flight,  relying  on  rotary  motion  of  a  complex  mechanical  system,  is 
predisposed  to  vibration.  While  the  helicopter  has  many  sources  of  vibrations 
the  main  rotor  system  generates  by  far  the  largest  magnitude  vibrations,  and  can 
render  the  vehicle  inoperable  if  left  unaddressed.  Thus,  proper  maintenance  to  reduce 
vibrations  is  essential  to  the  safe  operation  of  any  helicopter.  This  maintenance, 
however,  is  costly  and  time  consuming. 

Improving  the  maintenance  procedure  for  balancing  the  main  rotor  system  has 
been  an  area  of  active  interest  since  the  inception  of  the  helicopter.  However,  the  state 
of  the  art  in  rotor  balancing  still  requires  several  iterations  of  rotor  adjustments,  each 
necessitating  a  separate  test  flight  and  then  time  consuming  maintenance,  to  reduce 
the  vibrational  level  to  an  acceptable  amount.  The  intent  of  this  research  is  to  provide 
an  improved  rotor  vibrational  adjustment  methodology  that  significantly  reduces  the 
number  of  adjustment  iterations  required  to  reduce  main  rotor  vibrations. 

1.1  Present  Rotor  Vibrational  Reduction  Techniques 

Vibrations  of  the  largest  magnitude  in  the  main  rotor  system  are  primarily  the 
result  of  unbalanced  hub  mass  and  aerodynamic  forces.  These  forces  are  the  result 
of  the  individual  blades  exhibiting  unequal  lift  forces  as  they  perform  one  complete 
rotation  about  the  system  hub.  By  summing  the  resulting  forces  of  each  blade  about 
the  entire  azimuth  of  the  main  rotor  a  resultant  force  is  determined  that  nutates 
about  the  main  rotational  axis,  thus  creating  a  vibration.  Vibrational  reduction  is 
performed  by  adjusting  individual  main  rotor  blades  to  balance  out  the  lift  forces. 
Historically,  this  maintenance  was  referred  to  as  track  and  balance  as  the  general  idea 
was  to  adjust  the  rotor  blades  so  each  blade’s  tip  followed  the  same  path,  or  ’track’. 
With  each  tip  following  in  the  same  plane  of  rotation,  the  idea  was  that  each  blade 
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would  then  produce  the  same  lift.  This  is  based  on  the  assumption  that  if  the  blades 
are  identical  then  identical  loads  will  be  imparted  on  each  blade  that  follows  a  common 
track.  Identical  blades,  however,  do  not  truly  exist  and  thus  track  and  balance  is  not 
ideal  for  eliminating  vibrations.  Modern  rotorcraft  vibration  reduction  methods  rely 
on  adjustments  to  each  blade  to  reduce  measured  vibrations.  This  method  is  referred 
to  as  rotor  smoothing  as  blade  adjustments  are  made  to  ’smooth’  the  overall  system 
vibration  to  an  acceptable  level. 

1.2  Current  Advances  in  Rotor  Smoothing 

Rotor  smoothing  relies  on  a  defined  mapping  from  blade  adjustment  space  to 
system  vibration  space  in  order  to  reduce  overall  system  vibrations.  Historic  ro¬ 
tor  smoothing  methods  have  relied  on  empirically  determined  non-parametric  linear 
mappings  to  compute  a  rotor  balance  adjustment  solution  to  minimize  main  rotor 
vibrations.  Examples  of  such  an  approach  are  the  US  Army  Aviation  Vibration  An¬ 
alyzer  (AVA)  system  [1],  These  methods,  while  performing  better  than  simple  track 
and  balance,  produce  inaccurate  blade  adjustment  solutions  as  the  mappings  do  not 
completely  represent  the  system  under  test.  The  focus  of  current  research  has  been  to 
improve  the  system  response  mapping,  as  inaccurate  mappings  result  in  the  iterative 
adjustments  that  rotor  smoothing  is  plagued  with. 

One  suggested  approach  to  improving  the  system  mapping  is  to  replace  the  linear 
non-parametric  mapping  by  a  non-linear  neural  network  model.  This  model  has  been 
successfully  applied  to  the  US  Army  Vibration  Management  Enhancement  Program 
(VMEP)  program  [6]  .  While  this  method  has  been  shown  to  outperform  AVA,  it  relies 
on  a  non  parametric  mapping,  which  precludes  any  attempt  to  identify  modeling 
errors  by  reviewing  the  model  parameters  for  accuracy.  By  identifying  incorrectly 
identified  model  parameters,  the  VMEP  approach  has  the  ability  to  correct  inaccurate 
mappings.  This  will  ultimately  reduce  the  iterations  exhibited  by  current  approaches 
as  the  result  of  using  a  more  accurate  model  for  control  solution  development. 
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1.3  Research  Objectives 

It  is  the  intent  of  this  research  to  develop  a  linear  parametric  mapping  based 
method  for  rotor  vibrational  adjustment  in  which  individual  model  parameters  can 
be  verified  for  accuracy  and  likewise  adjusted  to  improve  model  accuracy.  To  achieve 
these  goals  several  steps  will  be  undertaken  as  enumerated  below. 

Objective  1:  Replace  the  non  parametric  mapping  of  the  rotor  system 
dynamics  with  a  parametric  approach  applicable  to  system  identification 
methods.  This  step  will  describe  the  methodology  defined  by  Wereley  [42]  to  model 
a  linear  system  with  time  periodic  coefficients  as  a  state  space  model,  in  a  manner 
similar  to  a  linear  time  invariant  system. 

Objective  2:  Adapt  the  Cramer-Rao  Bound  to  Validate  the  Parame¬ 
ters  of  the  Linear  Time  Parametric  Rotor  Model.  The  Cramer-Rao  Bound  is 
a  method  commonly  used  in  flight  testing  to  establish  the  accuracy  of  identified  pa¬ 
rameters  of  a  linearized  vehicle  model.  The  effort  of  this  step  is  to  develop  a  method 
to  adopt  the  Cramer-Rao  bound  to  the  parametric  model  defined  by  Objective  1. 

Objective  3:  Develop  a  Parametric  Main  Rotor  System  Model.  The 

next  step  in  achieving  the  proposed  research  goals  is  to  develop  a  parametric  time 
periodic  main  rotor  system  model  for  the  purposes  of  simulation.  The  model  will 
include  dynamically  actuated  and  fixed  pitch  linkages  for  each  blade  so  it  will  be 
possible  to  explore  both  static  and  dynamic  smoothing  cases. 

Objective  4:  Perform  System  Identification  of  the  Main  Rotor  Sys¬ 
tem.  This  step  in  the  proposed  research  goals  is  to  adapt  a  system  identification 
technique  to  determine  the  dynamics  of  the  main  rotor  in  both  hover  and  forward 
flight.  An  accurate  rotor  model  is  necessary  for  the  development  of  an  effective  vi¬ 
bration  controller,  as  will  be  done  as  the  final  objective  in  this  research. 

Objective  5:  Develop  a  Control  Methodology  for  Rotor  Vibration 
Smoothing.  The  final  objective  of  this  research  is  to  develop  a  robust  controller 
capable  of  attenuating  the  hub  vibrations  caused  by  aerodynamic  imbalances  of  the 
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rotor  system.  The  control  methodology  must  be  adaptable  to  the  periodic  nature  of 
the  helicopter  rotor  model  in  forward  flight.  This  will  use  the  verified  model  developed 
by  objectives  1-4. 

The  presentation  of  this  work  is  now  described.  The  next  chapter  presents  the 
historical  developments  of  rotorcraft  smoothing.  This  is  to  familiarize  the  reader  with 
the  successes,  but  more  importantly,  the  deficiencies  with  the  existing  rotor  vibration 
reduction  methods. 
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II.  Historical  Development  for  Research  in  Rotorcraft 

Smoothing 


The  practice  of  reducing  main  rotor  induced  vibrations  in  rotorcraft  has  been 
around  since  the  first  helicopters  were  developed.  This  chapter  presents  a  brief 
history  of  the  evolution  of  these  practices  with  the  intent  of  highlighting  their  suc¬ 
cesses  and  deficiencies.  Recommendations  to  overcome  the  existing  deficiencies  of  the 
methods  covered  in  this  chapter  are  presented  as  a  basis  for  the  work  this  research 
will  undertake. 

2. 1  Smoothing  of  Rotorcraft  Vibrations 

In  this  Section  we  discuss  the  effects  of  vibrations  on  a  rotorcraft  and  the  meth¬ 
ods  that  exist  to  reduce  them. 

2.1.1  The  Need  for  Rotorcraft  Vibration  Reduction. 

Helicopters,  as  with  any  rotating  system  generate  an  oscillatory  vibration  when¬ 
ever  the  forces  acting  on  the  system  are  imbalanced.  This  phenomenon  is  common¬ 
place,  as  all  who  have  driven  an  automobile  have  experienced  the  vibrations  of  an 
unbalanced  tire.  This  problem  is  not  only  an  annoying  disturbance  to  the  driver,  but 
if  left  untended  to  can  lead  to  costly  repairs.  In  the  case  of  the  automobile  the  proper 
maintenance  required  is  simply  to  rebalance  the  tire  by  determining  the  magnitude  of 
the  disturbance  force  and  at  what  phase  of  rotation  does  it  occur.  A  mass  that  gener¬ 
ates  an  equal  magnitude  force  is  placed  180  degrees  out  of  phase  with  the  disturbance. 
This  procedure  is  repeated  until  the  measured  vibrations  are  below  a  predetermined 
threshold. 

The  vibrations  exhibited  by  a  helicopter’s  rotating  blade  systems,  being  either 
the  main  or  tail  rotor,  are  similar  in  in  concept  to  the  automotive  example  above. 
These  vibrations  generate  the  highest  magnitude  vibrations  in  the  airframe  and  must 
be  balanced  out  to  prevent  crew  fatigue,  premature  airframe  fatigue,  and  catastrophic 
system  failures.  The  balancing  procedure  for  a  helicopter  is  in  principal  the  same  as 
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for  balancing  a  tire,  but  the  disturbance  forces  are  slightly  different  and  the  system 
complexity  is  greater  and  must  be  considered.  The  rotating  automotive  tire  has 
only  imbalanced  inertial  forces  to  cause  a  vibration,  whereas  the  helicopter  must  also 
contend  with  aerodynamic  force  imbalance  caused  by  the  individual  blades.  The 
aerodynamic  forces  can  also  vary  periodically  as  the  helicopter  transitions  from  a 
hover  state  to  a  forward  flight  state.  Aerodynamic  disturbance  forces  are  corrected 
by  adjusting  the  angle  of  attack  on  the  requisite  blades  as  required  to  meet  magnitude 
and  phase  requirements.  These  adjustments  can  be  made  by  changing  the  length  of 
a  pitch  linkage  or  by  adjusting  a  trim  tab  on  a  blade.  By  changing  the  length  of  the 
pitch  linkage  for  a  particular  blade  the  pitch  of  the  blade  is  adjusted  for  the  entire 
blade,  whereas  an  adjustment  to  a  trim  tab  adjusts  the  camber  of  the  blade  at  the 
portion  the  blade  in  which  it  is  attached.  This  procedure  is  generally  referred  to  as 
Track  and  Balance  or  more  commonly  Rotor  Smoothing. 

Rotational  vibrations  in  helicopters  generate  three  distinct  problems;  increased 
maintenance  downtime,  flight  crew  fatigue,  and  increased  vehicle  life  cycle  costs. 
Renzi  [33]  emphasizes  the  financial  impact  of  performing  maintenance  to  alleviate 
them.  He  points  out  that  the  cost  of  rotor  smoothing  maintenance  is  not  considered 
as  a  significant  cost  driver  during  vehicle  acquisition  but  becomes  one  of  the  most 
costly  operations  in  the  vehicle’s  lifespan.  A  more  efficient  method  of  rotor  smooth¬ 
ing  will  greatly  improve  the  operational  availability  of  the  vehicle,  crew  alertness,  and 
significantly  reduce  the  maintenance  cost  of  the  vehicle  over  it’s  lifetime. 

2.1.2  Rotor  Smoothing,  Track  and  Balance. 

There  are  two  common  terms  used  to  define  the  process  of  alleviating  helicopter 
1  per  revolution  rotor  vibrations  and  will  now  be  addressed  to  alleviate  any  confusion. 
The  first,  track  and  balance ,  is  a  more  historical  term  used  when  reducing  rotor  vibra¬ 
tions.  The  term  arises  from  the  earliest  methods  of  rotor  vibration  reduction  when 
the  blades  of  the  helicopter  were  adjusted  so  their  individual  blade  tips  are  aligned 
in  the  same  plane  of  rotation.  This  approach  noted  the  difference  in  vertical  spacing 
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between  the  individual  blade  tips,  usually  by  marking  them  in  different  colors.  The 
individual  blade  spacings  were  noted,  then  the  necessary  adjustments  to  each  blade 
are  made  to  bring  the  blade  tips  into  alignment  in  the  same  plane,  or  as  commonly 


referred  to,  in  the  same  ’track’.  An  example  of  this  is  depicted  in  Figures  2.1  and 


2.2.  This  was  the  first  method  to  attempt  to  balance  the  aerodynamic  forces  of  each 
blade  [33,35]. 


Figure  2.1:  Rotor  Blade  Tracking  [1], 


Color  mark 


Figure  2.2:  Image  of  Rotor  Blade  Split  [4]. 


The  second  term  used  when  defining  rotor  vibration  alleviation  is  smoothing. 
This  term  more  accurately  defines  the  modern  process  of  vibration  alleviation  as  the 
requirement  to  force  the  blades  to  track  in  the  same  plane  of  rotation  is  relaxed.  This 
process  usually  relies  on  performing  individual  blade  adjustments  that  will  reduce  the 
overall  chassis  vibrations  [39] .  Taitel  notes  that  since  the  reduction  of  rotor  vibrations 
are  the  main  goal  of  rotor  smoothing,  the  tracking  of  each  blade  may  not  be  perfect, 
as  the  adjustment  to  reduce  vibrations  may  result  in  a  split  track.  Rotor  smoothing 
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views  the  requirement  of  perfect  tracking  as  aesthetic  and  not  necessary  for  reduced 
vibration  levels. 

2.1.3  Cause  of  Rotor  Vibrations.  An  unbalanced  cantilevered  spinning 

rotor  can  emit  both  vertical  and  lateral  vibrations.  This  is  the  case  for  both  the  main 
rotor  and  tail  rotor  of  a  helicopter.  Lateral  vibrations  are  due  to  mass  imbalances, 
such  as  the  individual  blades  of  the  rotor  system  having  unequal  masses.  Vertical 
vibrations  are  due  to  the  individual  blades  in  the  rotor  system  having  unequal  lift, 
thus  causing  a  nutating  lift  vector  about  the  rotational  axis  of  the  rotor  system. 

Rotor  systems  emit  vibrations  across  an  infinite  frequency  band.  The  largest 
magnitude  vibrations  are  those  occurring  at  the  system  fundamental  frequency,  which 
is  once  per  revolution  (1/Per)  [5].  The  current  methods  that  exist  for  rotor  smoothing 
can  only  reduce  vibrations  at  the  fundamental  rotor  frequency.  This  is  due  to  current 
helicopter  control  systems  inability  to  command  anything  but  a  cyclic  control  at  the 
fundamental  frequency  of  the  rotor  system.  Higher  multiples  of  the  fundamental  rotor 
frequency  are  also  noticeable  sources  of  vibration  but  have  vibrations  that  are  orders 
of  magnitude  lower  than  the  fundamental  frequency.  These  vibrations  are  primarily 
due  to  the  harmonic  forces  generated  in  forward  flight  ,  which  are  due  in  part  to  the 
flexible  nature  of  the  rotor  system.  The  research  area  of  Higher  Harmonic  Control 
(HHC)  [12, 20]  is  addressing  this  problem. 

2.2  Rotor  Smoothing  Methods 

Rotor  smoothing  generally  is  the  process  of  adjusting  rotor  blade  properties  to 
reduce  the  vibrations  due  to  unbalanced  loads  across  the  rotor  disk.  This  Section  will 
review  the  adjustments  to  the  rotor  system  used  in  this  process. 

2.2.1  Rotor  Adjustment  Options.  Vibrations,  as  stated  previously,  are  due 
to  asymmetrical  forces  acting  on  a  spinning  rotor  system.  In  order  to  perform  any 
adjustments  to  correct  the  vibrations  the  control  inputs  that  are  available  to  a  heli- 


copter  must  first  be  know.  This  Section  will  cover  commonly  available  control  inputs 
on  modern  rotorcraft. 

2. 2. 1.1  Mass  Adjustment.  Mass  adjustments  are  necessary  to  correct 
the  lateral  vibrations  of  a  spinning  rotor.  The  rotor  system  provides  for  this  correction 
by  applying  masses  to  the  root  of  an  individual  blade  of  the  rotor  system,  as  seen  in 
Figure  2.3.  These  masses  are  usually  in  the  form  of  plates  or  pellets.  It  is  important  to 
note  the  effect  of  adding  mass  to  a  blade  can  be  replicated  by  removing  an  equal  mass 
from  the  opposite  blade,  or  opposite  blades  if  the  rotor  system  has  an  odd  number 
of  blades.  Adjusting  the  masses  on  rotor  blades  has  no  aerodynamic  effect  and  thus 
provides  no  direct  input  to  vertical  vibrations. 


Figure  2.3:  Blade  Weights  Applied  at  the  Blade  Root  [23]. 


9 


2. 2. 1.2  Pitch  Link  Adjustment.  In  order  to  affect  the  rotor’s  vertical 
vibrations  an  aerodynamic  input  must  be  imparted  on  the  system.  One  method  for 
this  is  by  adjusting  the  pitch  linkage  of  an  individual  blade.  A  pitch  linkage  is  a  rod  of 
adjustable  length  which  connects  the  blade  to  the  swashplate  of  the  helicopter.  This 
adjustment  is  seen  in  Figure  2.4. 


Figure  2.4:  Rotorcraft  Pitch  Linkage  Adjustment  [23]. 


The  pitch  linkage  controls  the  pitch  of  the  blade  it  is  connected  to.  Increasing  or 
decreasing  the  length  of  a  pitch  linkage  will  likewise  increase  or  decrease  the  angle  of 
attack  of  the  modified  blade.  The  angle  of  attack  of  the  blade  will  determine  the  lift 
the  blade  will  produce.  Thus,  for  positive  vertical  forces  to  be  reduced  on  the  rotor,  a 
blade  that  is  180  degrees  out  of  phase  of  the  disturbance  has  the  pitch  linkage  set  to 
increase  positive  pitch  to  impart  an  equal  and  opposite  force.  Alternately,  negative 
pitch  can  be  applied  to  a  blade  that  is  directly  in  phase  to  the  disturbance  to  achieve 
the  same  effect.  Pitch  linkages  on  modern  helicopters  are  not  dynamically  adjustable, 
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thus  any  change  in  pitch  length  must  occur  when  the  helicopter’s  rotor  system  is  not 
spinning.  Research  is  underway  to  provide  for  dynamically  adjustable  pitch  linkages. 
The  concept  is  seen  in  Figure  2.5. 


Figure  2.5:  Rotorcraft  Dynamic  Pitch  Linkage  [20]. 


Mannchen  [20]  and  Hwang  [12]  review  the  various  aspects  of  implementing  HHC 
via  dynamically  adjustable  control  linkages.  A  BO  105  has  been  modified  to  replace 
the  static  pitch  linkages  with  dynamically  adjustable  linkages,  as  seen  in  Figures  2.6 
and  2.7. 


2.2. 1.3  Trim  Tab  Adjustment.  Trim  tab  adjustments  are  another 
method  of  adjusting  the  rotor  system’s  vertical  vibrations.  A  trim  tab  is  a  flap-like 
device  attached  to  the  trailing  edge  of  a  rotor  blade.  An  example  of  a  trim  tab  is  seen 
in  Figure  2.8.  It  operates  in  much  the  same  manner  as  an  aileron  on  a  fixed  wing 
aircraft  by  changing  the  overall  camber  of  the  wing  section  to  which  it  is  attached. 
Thus,  reduced  camber  will  reduce  or  change  the  direction  of  lift  on  the  wing  at  the 
section  to  which  it  is  connected.  Likewise,  an  increase  in  camber  will  cause  an  increase 
in  lift  at  the  airfoil  section.  Trim  tabs  are  not  dynamically  adjustable,  but  similar 
research  as  discussed  by  Hwang  and  Mannchen  [12, 20]  allows  for  such  capability. 
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Figure  2.6:  BO  105  Rotor  Head  With  Dynamic  Pitch  Linkages  [20]. 

2.2.2  Historical  Methods.  Rotor  smoothing  methods  have  been  around  as 
long  as  helicopters  have  been  vibrating,  which  is  to  say  since  their  inception.  This 
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Figure  2.7:  Individual  Dynamic  Pitch  Linkage  on  a  BO  105  [20]. 
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Figure  2.8:  Rotor  Blade  Trim  Tab  at  Trailing  Edge  [33]. 


Section  will  briefly  cover  the  methods  used  from  the  1930’s  to  the  1980’s,  as  outlined 
by  Johnson  [15]. 


2.2.2. 1  Touch  Flag  Method.  The  touch  flag  method  was  the  forerunner 
in  rotor  track  and  balance  techniques.  This  method  was  devised  to  align  the  blade 
tips  in  the  same  plane  of  rotation,  or  ’track’.  This  method  employed  a  pole  mounted 
flag  with  which  the  operator  would  allow  to  come  in  contact  with  the  rotating  blade 
tips.  As  each  blade  tip,  which  was  colored  with  chalk  or  crayon,  struck  the  flag  the 
blade  gap  separation  was  then  noted,  thus  providing  track  information. 


2. 2. 2. 2  1960’s  Electro- Optical  Track  Adjustment.  In  the  1960’s  a 
method  to  measure  the  blade  track  using  electro-optical  equipment  was  devised  by 
Chicago  Aerial.  This  method  provided  for  very  accurate  measurement  of  blade  track 
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via  a  vertical  height  measurement  form  the  sensor  datum.  This  method,  however, 
only  provided  track  information  while  the  vehicle  was  on  the  ground.  A  second  prod¬ 
uct  was  devised  to  provide  in-flight  tracking  but  produce  lower  accuracy  results  as 
compared  to  the  the  high  accuracy  method. 

2. 2. 2. 3  1970’s  Strobe  Light  Tracking.  A  method  to  identify  blade 
track  information  while  the  vehicle  was  in  flight  was  devised  in  the  1970’s  using  a 
synchronized  strobe  and  reflective  targets.  The  targets  were  affixed  to  the  individual 
blade  tips.  This  method  required  the  operator  to  remember  individual  blade  positions, 
which  required  not  only  training  but  high  proficiency.  While  this  method  was  capable 
of  blade  tracking  while  the  vehicle  was  in  flight,  it  was  deemed  too  unreliable  due  to 
the  high  operator  proficiency  requirement  and  thus  abandoned. 

Rotor  smoothing  during  1970’s  moved  from  the  purely  static  mass  balance 
method  used  up  until  this  point  to  a  method  which  measured  rotational  vibrations. 
This  method  was  adopted  from  spin  balancing  techniques  for  large  industrial  blowers. 
The  data  from  the  accelerometers  was  synchronized  with  a  strobe  light  to  establish 
the  phase  relationship  of  the  vibration.  This  provided  the  operator  information  of 
how  much  and  where  to  place  mass.  This  method  did  not  address  blade  adjustments 
for  vibration  reduction. 

2. 2. 2. 4  1980’s  to  late  1990’s.  Rotor  vibrational  reduction  methods 
during  the  1980’s  began  to  adopt  a  mathematical  model-based  approach  by  using  a 
linearized  model  of  the  rotor  dynamics  to  determine  the  appropriate  rotor  adjust¬ 
ments.  This  approach  was  intended  to  eliminate  the  multi-step  process  of  previous 
rotor  track  and  smoothing  processes  by  identifying  the  relationship  between  rotor  ad¬ 
justments  and  vibrations  via  a  mathematical  model.  This  procedure  worked  well  as 
long  as  the  vehicle  is  a  close  match  to  the  mathematic  model  used  by  the  smoothing 
algorithm.  The  mathematical  models  used  primarily  non-parametric  transfer  function 
models.  As  the  mathematical  model  did  not  adjust  to  match  the  characteristics  of  the 
test  aircraft,  the  reliability  of  the  recommended  adjustments  were  valid  only  50  to  75 
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percent  of  the  time  [15].  This  type  of  algorithm  is  referred  to  as  a  non-learning  type 
as  the  model  is  rigidly  fixed  to  one  aircraft  configuration.  This  has  led  to  the  develop¬ 
ment  of  learning  algorithms  so  that  the  internal  model  can  adjust  to  match  the  test 
aircraft,  as  discussed  in  reference  [41]  and  in  Section  2. 2. 3. 2.  The  research  covered  in 
this  dissertation  primarily  focuses  on  the  improvement  of  learning  algorithms. 

2.2.3  Current  Methods.  Modern  developments  in  rotor  vibration  reduction 
are  focused  on  providing  an  adaptive  mathematical  model  which  can  ’learn’  the  indi¬ 
vidual  characteristics  of  the  the  vehicle  under  test.  These  methods  attempt  to  match 
the  input-output  relationship  of  test  data  to  that  of  the  mathematical  model  used  by 
the  smoothing  algorithm.  These  algorithms  are  based  on  both  linear  and  nonlinear 
models.  In  this  Section  we  will  briefly  describe  the  linear  and  nonlinear  algorithms 
used  to  smooth  rotor  vibrations.  The  cases  considered  here  are  for  a  single  main  ro¬ 
tor,  however  the  same  methodology  applies  for  the  tail  rotor  or  any  rotor  in  a  counter 
rotating  main  rotor  configuration,  such  as  a  Kamov  Ka-26. 

2.2.3. 1  Linear  Smoothing  Algorithms.  Linear  algorithms  were  the 
first  of  the  numerically  based  approaches  to  rotor  vibration  smoothing  to  appear  in 
the  helicopter  community.  This  approach  is  based  on  a  non-parametric  input-output 
relationship  of  the  change  in  vibrations  to  the  change  in  a  rotor  adjustment,  such 
as  changing  the  length  of  a  pitch  link  or  changing  the  angle  of  a  trim  tab.  The 
input-output  relationship  is  considered  linear  in  this  approach  and  thus  results  in  a 
simplified  transfer  function  model.  The  transfer  function  is  expressed  simply  as  a 
Least-Squares  Equation  2.1  [33] 


(x%]  =  AjVibration/  A  Adjustment^  (2-1) 

where  a%3  is  the  sensitivity  coefficient  matrix,  AjVibration  is  the  change  in  vibrations 
for  time  j  and  A  Adjustment  Si  is  the  change  in  adjustment  i  to  the  rotor.  It  is 
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important  to  note  that  the  linear  sensitivity  coefficients  are  determined  a  priori  to  any 
adjustment  calculations  through  a  series  of  flight  tests  at  select  altitudes  and  velocities 
within  the  helicopter’s  flight  envelope.  The  individual  elements  of  the  sensitivity 
coefficient  matrix  a  are  determined  by  making  incremental  changes  to  A Adjustinenti 
in  Equation  2.1.  This  process  is  repeated  for  all  incremental  changes  of  adjustments 
at  all  the  selected  points  within  the  helicopter’s  flight  envelope.  In  most  cases  in  linear 
smoothing  algorithms  the  sensitivity  coefficient  remains  unchanged  after  the  initial 
matrix  development.  This  approach  is  referred  to  as  a  non-learning  approach  as  the 
algorithm  does  not  adjust  the  sensitivity  coefficients  after  they  have  been  developed. 

Non-learning  algorithms  have  problems  with  accuracy  in  their  adjustment  cal¬ 
culations  due  to  error  in  the  individual  elements  within  the  sensitivity  coefficient 
matrix.  These  errors  arise  from  the  method  in  which  the  linear  coefficients  are  cal¬ 
culated,  namely  that  they  are  calculated  from  one  helicopter.  Individual  helicopters 
have  variances  in  their  input-output  relationships,  thus  the  sensitivity  coefficient  ma¬ 
trix  does  not  accurately  apply  to  all  helicopters.  Additionally,  error  arises  in  the 
sensitivity  coefficients  due  to  poor  reproducibility  of  the  vibration  data  at  each  flight 
condition.  An  individual  linear  coefficient  within  the  sensitivity  matrix  is  formed 
from  the  average  value  of  several  tests  at  the  same  flight  condition  and  same  adjust¬ 
ment  setting.  Wroblewski  [5]  indicates  that  the  statistical  variance  between  tests  at 
one  point  is  30%.  This  is  due  to  measurment  noise,  aircraft  gross  weight,  and  the 
interaction  of  aircraft  modes  that  are  weakly  dependent  on  rotor  adjustments  [5]. 

The  required  adjustment  is  determined  by  adjusting  Equation  2.1  as  seen  in 
Equation  2.2 


A  Adjustmenti  =  /  AjVibration  / .  (2.2) 

Once  again,  a  Least-Squares  approach  is  used  to  calculate  the  required  rotor  adjust¬ 
ment  necessary  to  minimize  the  measured  vibration. 


17 


2. 2. 3. 2  Artificial  Neural  Network  Smoothing  Algorithms.  Artificial 
neural  networks  are  another  method  of  describing  the  input-output  relationship  of 
a  dynamic  system.  Based  on  biological  neural  networks,  artificial  neural  networks 
mimic  the  functions  of  the  axion,  the  cell  body,  and  the  dendrites  as  seen  in  Figure 
2.9.  The  dendrites  act  to  accept  inputs  to  the  system,  the  cell  body  assigns  weightings 
to  inputs  and  then  passes  this  value  through  an  internal  transfer  function,  and  these 
signals  are  then  carried  out  of  the  cell  body  by  the  axion.  The  synapse  acts  as  the 
output  of  the  system  [21]. 


Figure  2.9:  Biological  Neural  Network  [21]. 

Artificial  neural  networks  allow  for  both  linear  and  nonlinear  transfer  functions 
within  the  ’’cell  body,”  thus  allowing  for  both  nonlinear  or  linear  mappings  of  input- 
to-output.  A  depiction  of  an  artificial  network  is  seen  in  Figure  2.10.  This  has  a 
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great  advantage  over  linear  transfer  function  models,  as  they  are  constrained  to  only 
linear  mappings  from  the  input  space  to  the  output  space.  Neural  networks  are  non- 
parametric  in  that  no  explicit  relationship,  such  as  a  multivariable  dynamic  model, 
is  needed  to  describe  the  relationship  of  helicopter  rotor  adjustments  to  vibrational 
output  [6].  The  model  is  therefore  represented  by  the  input-to-output  interactions  of 
the  network. 


Figure  2.10:  Artificial  Neural  Network  [21], 

Artificial  neural  networks  establish  the  input-to-output  mapping  through  a 
‘learning’  process  in  which  true  input  and  output  data  is  observed  and  then  repro¬ 
duced.  The  data  is  reproduced  by  first  passing  the  input  data  through  the  network 
and  then  comparing  the  network  calculated  output  to  the  true  output.  A  gradient  of 
the  magnitude  of  the  error  in  the  output  is  then  used  to  adjust  the  weightings  in  the 
individual  neurons  within  the  network  which  will  allow  the  network  to  then  calculate 
a  more  precise  output.  This  process  is  repeated  until  the  error  between  the  output  of 
the  neural  network  and  the  real  data  falls  below  a  predetermined  threshold. 

In  the  case  of  using  artificial  neural  networks  for  smoothing  main  rotor  vibrations 
a  set  of  training  data  must  first  exist.  This  set  of  data  is  referred  to  as  the  training 
set.  The  training  set  can  come  from  either  real  aircraft  or  simulation  of  a  linear  or 
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nonlinear  parametric  model.  In  the  establishment  of  the  training  set  a  range  of  rotor 
adjustments  is  performed  and  the  resulting  change  in  helicopter  vibration  is  then 
recorded  for  each  change.  Wroblewsky  [6]  indicates  the  training  set  must  be  selected 
to  exhaustively  cover  the  entire  flight  envelope  of  the  vehicle.  As  the  data  required 
to  adequately  train  a  neural  network  for  this  task  is  quite  large,  simulation  data  is 
typically  used  as  a  reasonable  first  approximation  [6].  As  more  flight  data  becomes 
available  the  training  set  is  modified  by  replacing  the  simulated  data  with  actual  flight 
data.  A  typical  training  set  will  require  between  20  to  30  flights. 

2.3  Smoothing  Performance  Of  Current  Methods 

Smoothing  performance  can  be  considered  in  terms  of: 

1.  Quality  of  final  rotor  adjustments 

2.  Number  of  rotor  adjustment  iterations  required  to  achieve  minimum  vibrations 

3.  Number  of  individual  rotor  adjustments  required  per  iteration 

4.  Amount  of  data  necessary  for  initial  non-parametric  mapping 

The  linear  non-parametric  mapping  method  presented  in  Section  2. 2. 3.1  and 
the  neural  network  based  method  presented  in  Section  2. 2. 3. 2  both  have  been  shown 
to  produce  adequate  adjustment  solutions,  thus  resulting  in  minimum  rotor  vibration 
[8,23,28-30].  In  this  comparison  there  is  no  clear  advantage  to  either  method  as  they 
both  arrive  at  an  acceptable  vibration  measurement  at  the  final  adjustment  solution. 
In  terms  of  speed  of  maintenance  there  also  seems  to  be  no  clear  leader.  Both  methods 
generally  require  two  to  three  smoothing  iterations  including  adjustments  and  test 
flights  to  reach  an  acceptable  level  of  rotor  vibrations.  This  is  interesting  as  the  neural 
network  method  was  reported  to  provide  a  solution  in  fewer  iterations  by  including 
the  system  nonlinearities  in  the  adjustment  to  vibration  mapping  [6].  Miller  [23] 
disputes  this  claim  by  demonstrating  that  rotor  adjustments  show  a  linear  mapping 
to  the  vibration  changes.  Furthermore,  the  number  of  changes  required  per  iteration 
appears  to  be  independent  of  method.  Both  methods  have  minimum  adjustment 
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optimized  solutions  available  and  perform  generally  the  same  [6,8].  Finally,  both 
methods  require  a  large  database  of  system  adjustments  and  resulting  vibrational 
data  to  build  the  initial  input-output  mapping.  In  this  area  there  appears  to  be 
no  clear  leader,  as  both  methods  can  use  either  flight  data  or  data  from  simulation 
to  create  the  initial  mapping.  Furthermore,  once  the  initial  mappings  have  been 
established  each  method  provides  no  method  to  adjust  them  to  an  individual  aircraft 
that  has  slightly  different  dynamics  due  to  wear,  environment,  and  variances  in  vehicle 
components.  Finally,  this  method  cannot  provide  adjustment  solutions  outside  of  the 
training  set. 

2.f  Current  Method  Deficiencies 

The  rotor  smoothing  methods  covered  up  to  this  point  represent  the  evolution 
of  the  state  of  the  art  in  this  area.  The  methods  have  emerged  from  flags  and  grease 
pencils  used  to  make  rotor  track  adjustment  methods  to  using  linear  or  nonlinear 
non-parametric  transfer  functions  used  to  compute  rotor  changes.  While  the  non- 
parametric  transfer  function  methods  produce  an  adequate  rotor  adjustment  solution 
to  reduce  vibrations  they  have  several  deficiencies  that  need  further  attention.  These 
deficiencies  are  listed  below: 

1.  Extensive  data  collection  required  before  any  use  of  system 

2.  Inability  to  adjust  mapping  after  establishment 

3.  Both  methods  are  non-parametric,  thus  it  is  impossible  to  check  the  individual 
parameters  for  accuracy 

This  Section  will  address  the  above  deficiencies  and  briefly  suggest  solutions  that  this 
work  will  be  based  on. 

The  first  deficiency  in  the  current  smoothing  methods  is  the  extensive  data 
collection  necessary  to  create  the  input-output  mappings  of  rotor  adjustments  to 
changes  in  vibration.  This  requirement  is  necessary  regardless  of  method  used,  as  the 
neural  networks  require  this  data  for  the  training  set.  As  stated  in  2. 2. 3. 2,  the  data 
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collection  can  require  up  to  30  test  flights,  depending  on  the  number  of  test  points 
used  to  create  the  mapping.  To  improve  the  state  of  the  art  in  rotor  smoothing 
methods  this  requirement  must  be  greatly  reduced. 

The  inability  to  adjust  the  input-output  mapping  after  it  is  created  is  the  second 
area  of  improvement  in  rotor  smoothing  algorithms.  Generally,  only  one  helicopter 
with  one  fixed  configuration  is  used  to  collect  the  initial  data,  which  is  then  applied 
to  the  entire  fleet.  Individual  aircraft  modifications  and  repairs  can  cause  differences 
in  the  input-output  mappings  as  the  vehicle  dynamics  are  different  from  the  vehicle 
used  to  create  the  original  mapping.  Thus,  the  accurate  reduction  in  vibrations 
requires  adjustments  based  on  the  current  vehicle,  not  the  one  from  which  the  mapping 
was  created.  The  second  improvement  to  rotor  smoothing  algorithms  is  to  allow  for 
changes  in  the  input-output  mapping  to  match  a  helicopter  undergoing  smoothing 
operations. 

It  is  difficult  to  verify  the  accuracy  of  Non  parametric  derived  input-output 
relationships  to  known,  or  truth  models,  as  they  do  not  contain  discernible  parameters 
such  as  lift  coefficients,  etc.  Thus,  after  the  initial  data  is  collected  and  the  mapping 
is  created  there  is  no  method  to  verify  the  validity  of  the  input-output  relationships 
except  further  empirical  testing  and  statistical  analysis.  A  third  improvement  to  rotor 
smoothing  is  to  incorporate  parametric  models  that  use  specific  parameters  that  can 
be  checked  for  validity. 

2. 5  Summary 

In  this  chapter,  we  exhibited  an  overview  of  previous  work  from  the  pertinent 
research  areas  and  noted  the  research  objectives  of  this  present  study  as  they  arose 
from  previous  works.  In  the  next  chapter,  the  objectives  for  the  proposed  research 
are  discussed  in  sequential  order. 
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III.  Scope  of  Research 


Rotor  smoothing  has  an  extensive  history,  as  seen  in  the  historic  view  taken  by 
the  previous  chapter.  Focusing  on  future  developments  in  rotor  smoothing,  this 
chapter  will  briefly  glimpse  at  the  research  objectives  of  this  work  and  outline  their 
significance  to  improvements  in  this  field. 

The  objectives  described  in  this  chapter  layout  the  framework  that  this  research 
will  use  to  produce  a  parametric  model-based  rotor  smoothing  algorithm.  In  principle, 
this  research  is  based  upon  a  proposed  rotor  smoothing  method,  which  works  as 
follows: 

1.  Perform  system  identification  to  populate  the  rotor  system  parametric  model. 

2.  Validate  the  accuracy  of  the  parameters  of  the  rotor  system  model. 

3.  Produce  a  vibration  control  solution  using  linear  optimal  methods. 

The  effort  of  this  research  is  to  redefine  how  rotor  smoothing  is  performed  by  evalu¬ 
ating  the  following  items: 

1.  Replace  the  non-parametric  mapping  of  the  rotor  system  dynamics  with  a  para¬ 
metric  approach  applicable  to  system  identification  methods.  This  will  be  done 
by  using  a  linear  time  periodic  system  modeling  approach  developed  by  Were- 
ley  [42]  that  produces  the  convenient  state  space  linear  operator  form. 

2.  Identify  a  method  that  is  capable  of  validating  parameters  of  the  parametric 
rotor  model  used  in  the  above  item. 

3.  Identify  a  System  Identification  methodology  that  is  compatible  with  a  rotor  in 
both  hover  and  forward  flight. 

4.  Produce  a  control  solution  that  reduces  rotor  vibration  levels  to  an  acceptable 
level  base  on  the  above  items. 

The  remainder  of  this  chapter  will  discuss  each  of  the  items  listed  above  as  objectives 
within  the  scope  of  this  work. 
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3.1  Scope  of  current  research 

The  primary  objective  of  this  research  is  to  develop  a  robust  rotor  smoothing  al¬ 
gorithm  based  on  frequency  domain  system  identification  methodology.  As  mentioned 
above,  this  method  will  be  compatible  with  a  helicopter  rotor  in  forward  flight,  which 
is  best  represented  as  a  linear  time  periodic  system  [12,42],  The  proposed  research 
outlined  below  is  to  extend  the  works  of  [12, 42]  by  introducing  a  theory  to  validate 
system  parameters  of  the  identified  model  based  on  the  Cramer- Rao  bound  [34] .  Ad¬ 
ditionally,  this  work  intends  to  address  the  feasibility  of  optimal  control  of  linear  time 
periodic  systems  for  vibration  smoothing  of  current  and  future  helicopter  main  rotor 
systems.  The  following  Sections  review  the  objectives  of  this  research  in  a  step  by 
step  approach. 

3.1.1  Objective  1:  Model  linear  systems  with  time  periodic  coefficients  as  a 
state  space  model.  This  Section  will  describe  the  methodology  defined  by  Wereley 
[42]  to  model  a  linear  system  with  time  periodic  coefficients  as  a  state  space  model,  in 
a  manner  similar  to  a  linear  time  invariant  system.  The  rationale  is  described  below. 

The  first  step  in  achieving  the  proposed  research  goals  is  to  develop  a  parametric 
main  rotor  system  model  for  the  purposes  of  simulation.  While  the  rotor  model  is 
inherently  nonlinear,  based  on  the  analysis  by  [23] ,  a  linear  model  is  deemed  adequate 
for  investigations  into  rotor  smoothing  and  thus  will  be  used  here.  This  is  advanta¬ 
geous  such  that  the  linear  model  approach  will  allow  the  use  of  developed  system 
identification,  analysis,  and  control  methodologies  for  linear  models.  Before  this  can 
begin,  however,  a  distinction  must  be  made  concerning  the  periodicity  of  the  system 
dynamics  and  how  to  address  them. 

The  rotor  system  model  proposed  above,  while  being  linear,  will  retain  the  time 
periodic  terms  in  the  system  dynamics,  control  input,  and  system  output  matrices. 
This  is  critical  in  order  to  retain  the  system  response  fidelity  in  forward  flight,  where 
the  time  periodic  terms  play  a  major  part  in  the  total  system  response.  Historically, 
these  terms  have  been  averaged  over  the  system  time  interval  T,  as  in  [20,31],  but  this 
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results  in  an  inaccurate  dynamics  model  and  ultimately  to  poor  control  development. 
This  averaging  is  due  to  the  inability  of  incorporating  the  time  periodic  terms  in 
the  state  space  form  of  the  linear  model.  To  explain  this  further,  it  is  important  to 
differentiate  between  a  model  with  Time  Invariant  dynamics  as  opposed  to  one  with 
Time  Periodic  dynamics. 

In  forward  flight,  a  helicopter  cannot  accurately  be  considered  as  a  Linear  Time 
Invariant  (LTI)  model  as  the  rotor  system  produces  periodic  aerodynamic  forces. 
This  requires  the  use  of  a  Linear  Time  Periodic  (LTP)  system  to  describe  the  rotor 
dynamics  in  forward  flight.  This  is  not  to  say  that  LTI  models  are  improper  to  use 
as  a  system  model  of  a  helicopter  rotor,  but  rather  that  the  flight  condition  must 
be  taken  into  account.  For  example,  Kvaternik  et  ah  [31]  successfully  used  a  LTI 
model  of  the  Bell-Boeiug  XV-15  tiltrotor  to  perform  real  time  system  identification 
and  control  to  reduce  rotor  vibrations  in  hover.  This  is  not  the  case  in  forward  flight, 
however,  as  Hwang  [12]  demonstrates  that  system  identification  based  on  LTI  models 
fails  due  to  poor  model  matching  as  compared  to  LTP  based  system  identification  of 
a  helicopter  in  forward  flight.  In  a  frequency  domain  sense,  this  is  due  to  the  inability 
of  the  LTI  system  to  account  for  the  sideband  power  generated  by  the  LTP  system,  as 
LTI  systems  can  only  account  for  power  at  the  excitation  frequency.  This  is  explained 
in  more  detail  in  Section  3. 1.1.1. 

3. 1.1.1  Frequency  Response  Differences  in  LTI  and  LTP  Systems.  As 
stated  above,  the  frequency  response  of  a  LTI  system  and  a  LTP  system  vary  in  output 
to  the  same  forced  input.  This  is  due  to  the  fact  that  a  LTP  system  can  change  the 
amplitude  and  phase  of  the  input  signal  in  addition  to  causing  frequency  translation 
of  the  input  signal.  LTI  systems  can  only  effect  the  amplitude  and  phase  of  the  input 
signal.  Hwang  states  [12]  the  frequency  translation  of  the  input  signal  in  LTP  systems 
is  represented  by  an  output  that  is  the  the  product  of  the  sinusoidal  excitation  signal 
and  a  Fourier  series  expansion  of  the  fundamental  frequency,  which  in  this  case  is  the 
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rotational  frequency  of  the  rotor  system.  An  example  of  this  is  shown  in  Figure  3.1. 
Linear  time  periodic  system  theory  will  be  covered  in  detail  in  chapter  V. 


Linear  Time-Constant  System: 

Sinusoidal  Input  Sinusoidal  Response 


Figure  3.1:  Multiharmonic  Response  of  an  LTP  System  [20]. 

Therefore,  a  LTP  system  will  produce  power  at  the  excitation  frequency  as  well 
as  sideband  power  due  to  the  fundamental  frequency  of  the  system.  It  is  clear  to  point 
out  that  any  system  identification  techniques  that  do  not  account  for  the  additional 
sideband  powers  will  result  in  inaccurate  results. 

3.1.2  Objective  2:  Adapt  the  Cramer-Rao  Bound  to  Validate  the  Parameters 
of  the  Linear  Time  Parametric  Rotor  Model.  This  Section  will  cover  the  neces¬ 
sity  to  develop  the  Cramer-Rao  Bound  to  validate  the  parameters  of  the  linear  time 
parametric  system  in  state  space  form,  as  in  the  case  of  a  helicopter  rotor  model. 
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System  models,  as  in  the  LTP  model  in  Section  3.1.1,  are  initially  described  by 
parameters  derived  from  mathematical  models  of  the  system  dynamics.  These  models, 
while  generally  good,  need  to  be  validated  so  that  they  accurately  describe  the  real 
system  upon  which  they  are  based  such  as  a  real  helicopter  rotor.  This  validation 
is  referred  to  as  system  identification ,  where  the  individual  system  parameters  are 
derived  from  data  collected  from  a  test  of  the  real  system,  such  as  a  flight  test. 

The  parameters  developed  by  system  identification  methods  define  a  linearized 
vehicle  model,  which  in  turn  defines  the  dynamic  characteristics  of  that  vehicle.  These 
models  are  in  turn  used  in  simulators,  control  system  design,  and  to  validate  wind 
tunnel  parameter  predictions.  Maine  and  Iliff  [34]  point  out  that  it  is  important  to 
remember  that  parameters  obtained  from  testing  are  only  estimates  and  not  exact 
values.  This,  unfortunately,  is  often  disregarded  and  rarely  are  the  parameters  ever 
verified  for  accuracy.  Maine  and  Lliff  [34]  further  state  that  if  accurate  parameter 
estimates  cannot  be  distinguished  from  worthless  estimates  all  estimates  must  be  as¬ 
sumed  to  be  of  questionable  accuracy.  It  is  for  these  reasons  that  parameter  validation 
methods  have  been  developed  for  LTI  systems  such  as  those  cited  in  [3,34] 

As  stated  in  Section  3.1  the  intent  of  this  research  is  to  develop  an  adaptive  con¬ 
trol  algorithm  to  alleviate  the  main  rotor  vibrations  for  the  purpose  of  rotor  smooth¬ 
ing.  The  measure  of  accuracy  of  the  model  parameters  of  the  identified  rotor  system 
are  of  critical  importance  in  this  application  as  they  will  determine  the  effectiveness  of 
the  implemented  vibration  controller.  To  implement  a  controller  based  on  an  unver¬ 
ified  model  may  have  disastrous  results  as  the  commanded  control  inputs  are  based 
on  a  model  that  does  not  match  the  true  rotor  system  dynamics. 

The  Cramer-Rao  bound  [34]  is  a  method  commonly  used  in  flight  testing  to 
establish  the  accuracy  of  identified  parameters  of  a  linearized  vehicle  model.  This 
measure  of  accuracy  is  based  on  the  Uncertainty  Ellipsoid  and  is  similar  to  other 
measures  of  accuracy  such  as  estimated  variance  and  standard  deviation.  The  Cramer- 
Rao  bound  is  the  same  as  these  methods  except  that  the  Cramer-Rao  bound  is  the 


27 


square  root  of  the  variance.  By  reviewing  the  Cramer-Rao  bound  for  each  parameter 
estimate  from  flight  test,  one  can  determine  the  data  accuracy  by  evaluating  the 
size  of  the  bounds  themselves.  The  Cramer-Rao  bound  established  the  standard 
deviation  of  an  identified  model  parameter,  and  therefore  large  bounds  indicate  a 
poor  estimation  performance.  This  method  provides  greater  insight  into  the  accuracy 
of  estimated  parameters  than  simply  reviewing  the  scatter  of  parameter  estimates 
from  flight  test  data  at  each  point,  as  seen  in  Figure  3.2.  Upon  review  of  Figure  3.2, 


O  Flight  data  solid  symbols  indicate  derivatives 

I  Cramer-Rao  bound  obtained  immediately  before 

-  Fairing  aircraft  departures 

Dashes  indicate  less  certain  fairings 


a,  deg 


Figure  3.2:  Example  of  Cramer-Rao  Bounds  for  Parameter  Estimates  [34], 


one  can  see  the  Cramer-Rao  bound  shows  that  data  below  -10  degrees  is  unreliable 
even  as  the  data  scatter  is  not  very  large.  Likewise,  at  +10  degrees  /alpha  the  Cramer- 
Rao  bound  matches  the  spread  of  the  collected  data.  In  this  case  the  parameter  is 
poor  based  upon  both  data  scatter  analysis  and  the  Cramer-Rao  bound.  This  shows 
that  the  Cramer-Rao  bound  can  detect  both  small  and  large  data  scatter  cases  where 
parameter  estimation  is  poor.  Thus,  a  review  of  the  Cramer-Rao  bounds  for  estimated 
vehicle  models  parameters  allows  for  an  evaluation  of  whether  the  parameters  have 
been  estimated  properly,  and  thus  whether  they  should  either  be  re-evaluated  or 
adjusted  by  the  engineer  before  they  are  used  in  the  vibration  control  model. 


This  research  will  develop  this  method  further  by  allowing  for  Cramer-Rao 
bounds  to  be  identified  for  LTP  model  parameters. 

3.1.3  Objective  3:  Develop  a  Parametric  Main  Rotor  System  Model.  The 
system  identification  an  control  processes  used  in  the  proposed  method  of  rotor 
smoothing  will  require  a  linear  time  periodic  model  of  a  helicopter  rotor.  This  Sec¬ 
tion  will  develop  a  linear  model  of  a  helicopter  rotor  system  which  incorporates  time 
periodic  system  coefficients  to  accurately  describe  the  system  in  forward  flight. 

The  model  will  include  dynamically  actuated  pitch  linkages  for  each  blade  so 
it  will  be  possible  to  explore  dynamic  smoothing  cases  allowed  by  using  an  optimal 
control  approach.  The  flap  equations  of  motion  for  a  rigid  blade  will  be  derived. 
Candidate  models  to  be  considered  are  those  based  on  the  work  presented  by  Johnson 
[15]  and  and  Webb  [37].  This  research  will  consider  main  helicopter  rotor  in  forward 
flight  and  in  hover,  however  the  emphasis  will  be  on  a  helicopter  in  forward  flight. 
The  inflow  will  be  simplified  to  consider  a  uniform  inflow  model. 

3.1.4  Objective  j:  Perform  System  Identification  of  the  Main  Rotor  System. 
The  next  step  in  the  proposed  research  goals  is  to  adapt  a  system  identification  tech¬ 
nique  to  determine  the  dynamics  of  the  main  rotor  in  both  hover  and  forward  flight. 
An  accurate  rotor  model  is  necessary  for  the  development  of  an  effective  vibration 
controller.  Therefore,  the  linear  time  periodic  model  proposed  in  Section  3.1.3  will  be 
used  for  this  step.  The  choice  of  system  identification  techniques  is  constrained  due 
to  the  inherent  dynamics  of  a  helicopter.  As  stated  in  Section  3.1.1,  in  forward  flight 
a  helicopter  cannot  accurately  be  considered  as  a  linear  time  invariant  model  as  the 
rotor  system  produces  periodic  aerodynamic  forces.  This  requires  the  use  of  a  linear 
time  periodic  system  to  describe  the  rotor  dynamics  in  forward  flight. 

Candidate  System  Identification  Methods.  System  identification 
methods  generally  fall  into  two  categories  for  linear  systems,  Frequency  Domain  meth¬ 
ods  and  Time  Domain  Methods.  Frequency  domain  methods  dominated  the  early 


29 


history  of  system  identification  with  non-parametric  transfer  function  analysis.  This 
method  was  improved  upon  by  the  developments  of  Cooley  and  Tukey  [13]  with  the 
Fast  Fourier  Transform(FFT),  which  greatly  decreased  the  testing  time  required  to 
generate  the  model.  Frequency  domain  methods  are  popular  in  system  identification 
as  they  have  the  ability  to  reject  low  frequency  drift  and  DC  bias  in  sampled  data. 
However,  Juang  [17]  states  that  frequency  domain  methods  have  lost  their  popularity 
in  control  model  identification.  Modern  control  methods  rely  on  parametric  models 
such  as  state-space  models.  Time  domain  methods  are  capable  of  generating  para¬ 
metric  models  and  for  this  reason  it  is  the  leading  system  identification  method  used 
for  the  development  of  control  model  identification.  Time  domain  methods  lack  some 
of  the  desirable  features  of  frequency  domain  methods  such  as  reduced  amounts  of 
sampled  data.  Additionally,  time  domain  methods  do  not  inherently  have  the  ability 
to  handle  bias  and  low  frequency  drift.  These  problems  can  be  addressed,  but  require 
the  use  of  additional  algorithms  and  thus  are  not  as  efficient  as  frequency  domain 
methods.  Recently,  several  researchers  [12, 42]  have  improved  frequency  domain  sys¬ 
tem  identification  methods  by  developing  a  method  to  identify  model  parameters  to 
match  the  frequency  response  characteristics.  Hwang  [12]  has  developed  a  parametric 
frequency  domain  identification  method  for  LTP  systems  which  is  directly  applica¬ 
ble  to  control  model  identification.  This  method  will  be  adapted  to  the  main  rotor 
smoothing  problem  of  this  research  due  to  its  inherent  ability  to  handle  LTP  methods 
and  generate  a  parametric  model. 

3.1.5  Objective  5:  Develop  an  LTP  Optimal  Control  Methodology  for  Rotor 
Vibration  Smoothing.  The  final  objective  of  this  research  is  to  develop  a  robust 
controller  capable  of  attenuating  the  hub  vibrations  caused  by  aerodynamic  imbal¬ 
ances  of  the  rotor  system.  The  general  idea  behind  attenuating  rotor  vibrations  is 
simply  inducing  equal  magnitude  forces  and  moments  into  the  rotor  system  that  are 
180°  opposite  in  phase  of  the  disturbance  force.  This  produces  destructive  interfer¬ 
ence  of  the  vibrations  caused  by  the  imbalance  loading  on  the  rotor  system,  effectively 
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eliminating  the  vibration  felt  at  the  rotor  hub  if  performed  precisely.  It  is  the  matter 
of  model  precision  that  the  previous  Section  3.1.3  dealt  with,  which  is  essential  to 
produce  a  linear  model  for  effective  controller  development.  The  research  of  Objective 
5  will  explore  the  use  of  active  control  to  arrive  at  a  steady  state  control  solution  to 
reduce  or  eliminate  rotor  induced  1/rev  vibrations.  The  following  Section  will  briefly 
discuss  candidate  control  methodologies  in  which  the  vibration  attenuation  system, 
or  simply  rotor  smoothing,  will  be  based  upon. 

3. 1.5.1  Candidate  Control  Methodologies.  Several  control  methodolo¬ 
gies  currently  exist  in  the  area  of  rotorcraft  control  including  Optimal  [7,36,40],  Modal 
Control  [37],  and  Periodic  system  control  [19,27].  Due  to  the  constraints  imposed  by 
the  modeling  requirements,  namely  that  the  model  must  be  represented  as  LTP  for 
accuracy,  only  methods  that  are  capable  of  handling  LTP  models  will  be  considered. 

Optimal  control  has  successfully  been  used  in  active  rotorcraft  vibration  control, 
as  sited  in  the  references  above.  Namely,  linear  optimal  methods  are  of  interest  to 
explore  as  they  inherently  use  parametric  models  for  the  controller  development.  This 
is  especially  attractive  as  the  intent  of  this  research  is  to  both  identify  and  verify  a 
parametric  linear  model.  The  application  optimal  LTP  control  has,  however,  only 
focused  on  active  vibration  reduction  for  an  already  smoothed  rotor  system.  This 
research  will  explore  the  usage  of  optimal  LTP  control  for  the  purpose  of  producing 
a  smoothed  rotor  system  by  producing  a  steady  state  control  solution.  As  this  is  the 
case,  only  the  1  per  rev  vibrations  will  be  considered  in  terms  of  a  disturbance  input. 
All  higher  multiples  of  disturbance  input  vibrations  will  be  disregarded  as  this  research 
is  not  focusing  on  higher  harmonic  control.  Steady  state  input  will  be  possible  for 
the  periodic  system  by  considering  all  control  inputs  in  the  fixed  coordinate  system 
of  the  rotor. 


This  chapter  presented  an  overview  of  the  overall  methodology  for  the  proposed 
research,  including  a  step-by-step  road  map  of  the  sequential  research  objectives.  The 
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next  chapter  will  introduce  the  mathematical  concepts  of  the  LTP  system,  which  will 
serve  as  the  basis  upon  this  work  will  be  built. 
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IV.  Mathematic  Foundations  of  Linear  Time  Periodic 

Systems 

inear  models  can  be  assumed  to  have  time  invariant  or  time  periodic  dynamics, 
as  discussed  in  Chapter  III.  While  the  principals  of  linear  time  invariant  models 
are  well  understood,  many  of  the  mathematical  foundations  of  linear  time  periodic 
systems  are  not  in  the  engineering  community.  Therefore,  this  Chapter  will  present 
several  of  the  required  mathematical  elements  of  linear  time  periodic  systems,  as  they 
will  be  used  extensively  in  later  Chapters  to  develop  the  LTP  system  model. 

4-1  The  Fourier  Series 

Harmonic  system  excitation,  in  the  case  of  input  signals  to  a  linear  transfer 
function,  repeat  at  a  time  interval, [0,  T],  The  term  T  is  referred  to  as  the  system 
period  and  is  defined  as 

_  2t r 
T  =  — 

UJp 

where  top  is  the  system  Fundamental  frequency  or  Pumping  frequency  (a ;p)  .  Any 
function  that  repeats  itself  over  the  time  interval  [0,  T]  is  known  as  a  periodic  function. 
Meirovitch  [22]  states  that  any  periodic  function  satisfies  the  relation  of  the  type 

f(t)  =  f(t  +  T)Vfe».  (4.2) 

By  the  definition  stated  in  Equation  4.2  it  is  evident  that  the  even  and  odd  trigono¬ 
metric  functions,  sinnt  and  cos  nt  ,  (n  =  1,2, . . .),  are  periodic  over  the  period  T 
as  sin  nt  =  sin  n{t  +  T)  and  likewise  for  cos  nt  when  up  —  1. 

An  example  of  a  periodic  system  excitation  is  a  sinusoidal  function  as  seen  in 
Figure  4.1.  From  this  graphic  it  is  clear  that  the  input  function  is  both  periodic  over 
T  and  purely  sinusoidal  over  that  period. 


(4.1) 
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It  is  important  to  note  that  a  harmonic  function  can  undergo  a  motion  that  is 
not  purely  sinusoidal  while  cycling  through  a  period  T.  The  saw  tooth  wave  function 
seen  in  figure  4.2  represents  an  example  of  a  non-sinusoidal  harmonic  function  that  is 
periodic  with  period  T.  Thus,  the  purely  sinusoidal  input  signal  is  a  general  case  of 
the  harmonic  input  signal.  Non-sinusoidal  functions,  while  still  periodic,  can  be  rep¬ 
resented  as  a  linear  combination  of  harmonic  functions.  This  combination  is  referred 
to  as  the  Fourier  Series  [22], 

To  develop  the  definition  of  the  Fourier  series  the  concept  of  linear  independence, 
orthogonality,  and  orthonormal  function  sets  must  first  be  defined.  These  subjects 
will  now  be  covered. 

4-1.1  Linear  Independence,  Orthogonality,  and  Orthonormal  Function  Sets. 
Over  the  interval  [0,  T]  a  set  of  functions  :  y  =  1,  2, ...  is  said  to  be  an  orthogonal 
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Figure  4.2:  Periodic  non-sinusoidal  input  signal  [22], 


set  of  functions  defined  on  [0,  T\  if 


ipr(t)ips(t)dt  =  0  V  r  7^  s  . 


(4.3) 


Furthermore,  the  set  of  functions  -0r(t)  :  r  =  1,2, ...  is  said  to  b eorthonormal  if 


il>r(t)dt  =  1  V  r  G  N 


(4.4) 


For  any  function  within  the  orthonormal  set,  the  Kronecker  Delta  (5)  can  be 
defined  as 


,ipr(t)'ips(t)dt  =  SrsW  r,s  . 


(4.5) 


Finally,  considering  the  set  of  functions  if)r(t)  :  r  —  1,  2, . . .  and  a  set  of  constants 
cr  G  3ft  :  r  =  1,2,...  where  cr  are  not  exclusively  null.  Then,  a  homogeneous  linear 
relationship  4.6 

n 

'EcrMt)  7^0  (4.6) 

7’=1 

denotes  the  function  set  ?/;r  is  defined  linearly  independent.  Likewise,  if  Equation  4.6 
is  equal  to  zero  then  the  set  ?/v  is  defined  as  linearly  dependent.  It  is  important  to  note 
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that  any  set  that  is  orthogonal  is  also  linearly  independent  by  noting  that  Equation 
4.3  by  definition  satisfies  the  properties  of  Equation  4.6. 


4-1.2  Trigonometric  and  Complex  Forms  of  the  Fourier  Series.  A  piecewise 
linear  function  f(t)  can  be  approximated  by  considering  a  orthonormal  set  f)r  in  the 
linear  relation  of  Equation  4.7. 

n 

■  (4.7) 

r=  1 

Any  function  set  {fr  :  r  =  1 . . .  n}  that  approximates  the  piecewise  continuous 
function  in  the  Equation  4.7  is  termed  complete  [22],  A  well  known  complete  function 
set  [22]  over  the  interval  [0,  27r]  is  seen  in  4.8. 


1  sint  cost  sin2t  cos2t  sinSt 

7^’  7^’ 


Using  the  complete  orthonormal  set  defined  in  Equation  4.8  every  continuous  function 
f(t)  defined  over  the  interval  [0,  27r]  can  be  represented  by  the  Fourier  Series ,  as  seen 
in  Equation  4.9. 


m 


00 

+  ^^(arcos  ru0t  +  brsinru0t ) 

r=  1 


(4.9) 


noting  oo0  is  referred  to  as  the  fundamental  frequency  .  This  form  of  the  Fourier 
series  is  referred  to  as  the  Trigonometric  form  of  the  Fourier  Series.  The  coefficients 
ar,  br  (r  =  1,2,.. .)  in  Equation  4.9  are  referred  to  as  Fourier  coefficients.  These 
coefficients  are  defined  as 


ar 

br 


—  /  f  (t)  cos  ru0t  dt,r 


-  /  f(t)  smru0t  dt,r 


0,1,2,... 
1,2,...  . 


(4.10) 
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The  Fourier  coefficients  define  the  contribution  of  the  r—1,2,. . .  trigonometric  terms 
cosruj0t  and  sinruj0t  to  the  approximation  of  the  piecewise  continuous  function. 

The  Fourier  series  can  also  be  expressed  as  an  exponential  series.  This  form, 
while  mapping  directly  to  the  trigonometric  form,  is  used  as  it  more  easily  expresses 
a  harmonic  excitation  in  exponential  as  compared  to  the  trigonometric  form.  The 
exponential  form  of  the  Fourier  Series  can  be  represented  as  follows 

OO 

/(f)  =  CreirUot  (4.11) 

r=— oo 

where  the  complex  Fourier  Coefficients  are  represented  as 

Cr  =  f(t)e~irUotdt,  r  =  0,±1,±2,...  .  (4.12) 

While  both  the  trigonometric  and  exponential  forms  of  the  Fourier  series  are  presented 
in  this  section,  this  work  henceforward  will  extensively  use  the  exponential  form. 

The  Fourier  series  in  Equation  4.9  represents  an  approximation  of  the  piecewise 
continuous  function  having  infinite  dimension.  In  most  cases  an  infinite  approximation 
is  not  feasible  and  must  therefore  truncate  the  number  of  terms  in  the  series.  The 
number  of  terms  included  in  the  Fourier  series  is  determined  by  the  level  of  accuracy 
required  to  reproduce  the  piecewise  continuous  function  it  is  approximating. 

4-2  Eigenvalues  and  Eigenvectors 

The  study  of  linear  systems  rely  heavily  on  the  matrix  exponential  function  ( eAt ) 
when  developing  solutions  to  differential  Equations.  These  solutions  can  be  used  to 
evaluate  the  system  characteristics  such  as  stability  and  frequency  response.  Key  to 
these  solutions  are  the  eigenvalues  and  eigenvectors  of  the  system.  This  section  will 
define  the  properties  of  these  terms  as  they  are  used  extensively  in  this  work. 
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Starting  with  the  definition  of  a  first  order  system  of  Equations  in  matrix  form 


du 

dt 


=  Au 


(4.13) 


where  the  constant  matrix  A  is  of  dimension  n  x  n.  Strang  [38]  notes  a  solution  to 
the  linear  Equation  4.13  as 


u{t)  =  ext(  .  (4.14) 

Now,  directly  substituting  Equation  4.14  into  4.13  produces  the  Equation 

A(  =  X(  .  (4.15) 

The  results  of  the  form  of  Equation  4.15  can  be  transformed  into  the  commonly 
recognised  form 


[XI-A](  =  0  (4.16) 

where  the  the  scalar  A  is  termed  the  Eigenvalue^ A)  and  the  vector  (  is  termed  the 
Eigenvector  (£)  of  the  matrix  A.  For  clarification  the  matrix  /  is  the  identity  matrix  of 
dimension  n.  Reid  [32]  notes  that  considering  the  scalar  first,  it  is  seen  that  A  must  be 
a  real  or  complex  value  such  that  the  coefficient  matrix  [A I  —  A]  is  singular,  thus  not 
invertible.  If  the  matrix  is  not  singular  then  the  trivial  solution  for  the  eigenvector 
would  be  the  zero  vector. 

Reid  [32]  further  notes  that  one  condition  for  the  matrix  [XI  —  A] ,  or  Resolvent 
Matrix  ([A/  —  A]),  to  be  non-singular  is  that  the  determinant  be  zero.  Thus,  the 
determinant  of  the  resolvent  matrix  yields  the  characteristic  polynomial  (A(X)),  an 
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nth  order  polynomial  in  A, 


A(A)  =  det[XI-A] 

—  A"  +  an_iA”  1  +  •  •  •  +  aiA  +  no  (4.17) 

where  a*  G  C.  The  eigenvalues  of  the  resolvent  matrix  are  the  distinct  roots  of  the 
characteristic  polynomial. 

4-3  Singular  Value  Decomposition 

One  method  of  determining  the  eigenvalues  and  eigenvectors  of  a  linear  system  is 
to  use  the  Singular  Value  Decomposition  (SVD).  Burl  [2]  defines  the  SVD  as  a  matrix 
factorization  that  provides  the  principal  gains  along  with  the  resultant  input  and 
output  direction  vectors  of  the  the  given  input-output  relationship  under  evaluation. 
The  SVD  is  essentially  a  similarity  transformation ,  which  is  defined  by  Wereley  [42] 
as  in  the  following. 

Definition  4.3.1  (Similarity  Transformation).  A  matrix  A  G  37nXn  can  be  reduced  to 
a  diagonal  matrix  A  by  the  similarity  transformation  A  =  WAV  if  and  only  if  A 
has  a  linear  independent  set  of  n  right  eigenvectors.  These  right  eigenvectors  are  the 
columns  of  V,  with  corresponding  left  eigenvectors  that  are  the  rows  of  W  =  V~x , 
and  corresponding  eigenvalues  that  are  the  diagonal  entries  of  A.  Also,  the  original 
matrix  A  can  be  reconstructed  from  its  eigenvectors  and  eigenvalues  as  A  =  VAW . 

Thus,  Burl  [2]  states  the  Singular  Value  Decomposition  is  defined  as  a  matrix 


A  =  VAW  (4.18) 

where  V  G  CnyXny  and  W  G  Cn“x"u  are  unitary  matrices,  noting  Cmxn(CmX") 
defines  a  set  of  complex  matrices.  Unitary  matrices  are  defined  as 
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Definition  4.3.2  (Unitary  Matrices).  A  matrix  V  having  the  property 


VfV  =  VVf  =  I  VU  E  Cnxn 


(4.19) 


is  thus  considered  a  unitary  matrix.  The  matrix  transpose  is  represented  by  the  symbol 

t- 

The  column  vectors  V  and  W  are  the  left  and  right  singular  vectors  of  A, 
respectively.  Furthermore,  the  matrix  A  E  5 lnvxnu  contains  the  eigenvalues  of  A,  or 
singular  values(cr),  as  seen  in  4.20 


A 


on 


0 


0 


(Jr) 


(4.20) 


where  the  dimension  ny  =  nu  =  np.  Burl  [2]  defines  similar  singular  value  matrices  of 
unequal  dimensions  y,u.  As  a  note  the  singular  values  appearing  on  the  diagonal  of 
A  appear  in  order  of  highest  to  lowest  value  such  that  oq  >  cr2  >  •  •  •  crp  >  0. 

By  performing  the  singular  value  decomposition,  Burl  [2]  states  that  one  can 
clearly  see  how  a  matrix  operates  on  a  vector,  or  input.  The  singular  values  represent 
the  range  of  matrix  gains,  which  are  referred  to  as  the  system  principal  gains.  They 
provide  a  range  of  maximum  and  minimum  gains  over  a  range  of  frequencies.  The 
right  singular  vectors  define  which  inputs  produce  the  maximum  and  minimum  gains, 
where  as  the  left  singular  vectors  define  which  outputs  result  when  the  maximum  and 
minimum  gains  are  achieved  [2]. 


4-4  The  Toeplitz  Transformation 

Many  of  the  modern  control  theories  have  been  based  upon  a  state  space  repre¬ 
sentation  of  a  linear  operator  form.  In  the  case  of  linear  operators  that  are  expanded 
by  a  Fourier  series  the  Toeplitz  transformation  is  used  when  converting  to  state  space 
form.  This  transformation  allows  the  algebraically  simplified  representation  of  the 
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state  space  matrix  representation  while  accounting  for  the  Fourier  series  components. 
This  section  will  define  the  Toeplitz  transform  and  the  properties  that  define  it,  as 
defined  by  Wereley  [42], 

Definition  4.4.1  (Toeplitz  Transform).  Let  the  T  periodic  matrix  Aft)  be  described 
by  the  absolutely  convergent  Fourier  series 

OO 

A(t)  =  J2  Anejnu)pt  (4-21) 

n= — oo 


where  Up  is  the  fundamental  frequency.  The  Toeplitz  Transform,  T  of  Aft),  TfAft)}, 
maps  the  set  of  complex  Fourier  coefficients  ,  {An\n  €  hi}  into  the  doubly  infinite 
block  Toeplitz  matrix,  A,  as  seen  in  Equation  f.22 


T{A}  =  A 


Ao 

A_! 

A—2 

1 

CO 

A- 4 

A1 

Ao 

A_! 

A— 2 

A_3 

A2 

A, 

Ao 

A_! 

A_  2 

A3 

A2 

A1 

Ao 

A_! 

a4 

A3 

A2 

A1 

Ao 

(4.22) 


The  Toeplitz  transform,  as  seen  above,  provides  a  convenient  method  of  trans¬ 
forming  periodic  differential  Equations  into  matrix  form.  Also  several  useful  proper¬ 
ties  defined  by  Wereley  [42]  associated  with  the  Toeplitz  transform  are  listed  below. 

Theorem  4.4.1  (Multiplicative  Property  of  the  Toeplitz  Transform).  The  Toeplitz 
transform  of  the  product  of  two  Toeplitz  matrices,  Aft )  and  Bft),  is  the  product  of  its 
transforms 


T{A  B}  =  T{A}T{B}  . 


(4.23) 
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Theorem  4.4.2  (Additive  Property  of  two  Toeplitz  Matrices).  The  Toeplitz  trans¬ 
form  of  the  sum  of  two  T  periodic  matrices  of  the  same  dimension,  Aft )  and  Bft),  is 
the  sum  of  its  transforms 

T{A  +  B}  =  T{A}  +  T{B}  .  (4.24) 

Definition  4.4.2  (Unitary  Matrices).  Consider  the  absolutely  convergent  series 

OO 

Aft)  =  X  jnuJpAnejnuJpt 

n=— oo 

as  A  G  (L2[0,  T],  Cnxn).  The  Toeplitz  transform  of  the  matrix  A  is  defined  as 

T{A}  =  A FT  {A}  -  T{A}  Af  (4.25) 

where  the  matrix  A f,  which  contains  multiples  of  the  pumping  frequency  up  is  defined 
as 


Af 


blkdiag{j7iujpI}\/  n  G  Z 
jujp  0 


Vn  G  Z  . 


0  j?iujp 


(4.26) 


This  Chapter  presented  an  overview  in  the  mathematical  preliminaries  necessary 
to  construct  the  foundation  of  linear  time  periodic  systems.  The  next  Chapter  will  de¬ 
velop  the  linear  time  periodic  system  theory  and  define  the  state  space  representation 
of  such  a  system. 
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V.  Linear  Time  Periodic  System  Theory 


This  Section  will  describe  the  formation  of  the  Linear  Time  Periodic  system  in 
terms  of  the  state  space.  The  linear  time  periodic  state  space  representation 
is  essential  to  later  chapters  in  system  parameter  validation  and  control.  In  the 
upcoming  chapters  linear  time  periodic  analogues  of  system  parameter  validation  and 
optimal  control  methods  will  build  upon  existing  state  space  based  LTI  methods. 
These  efforts  will  first  require  a  linear  time  periodic  state  space  operator,  which  this 
chapter  will  detail. 

5.1  The  Continuous  Time  Invariant  Linear  State  Model 

The  linear  continuous  time  state  model  is  a  matrix  representation  of  several 
first  order  linear  differential  Equations.  This  matrix  representation  presents  a  linear 
mapping  of  the  input  to  output  response  of  the  linearised  system  dynamics.  The 
standard  representation  of  the  state  space  model  is  seen  in  Equations  5.1  and  5.2. 


x{t)  =  A(t)x(t)  +  B(t)u(t) 

(5.1) 

y(t )  =  C(t)x(t)  +  D(t)u(t) 

(5.2) 

In  this  form  the  states  of  the  system  are  represented  as  the  vector  x(t)  G  Mnx  which  act 
upon  the  system  State  or  Plant  matrix,  A(t)  G  MTl:cXn:':.  The  vector  u(t)  G  Mn“  rep¬ 
resents  the  system  input  which  acts  upon  the  system  Control  matrix  B(t)  G  Mn-X'xn,i, 
as  seen  in  the  state  equation  5.1.  The  system  output,  y(t),  as  seen  in  the  output 
equation  5.2,  is  composed  of  the  the  state  vector  x(t)  acting  upon  the  system  Output 
matrix  C(t)  G  WlyXnx  and  control  vector  u(t)  acting  upon  the  Feed  Forward  matrix 
D(t)E  WlyXnu.  It  is  important  to  note  that  the  A(t),  B(t),  C(t)  and  D(t)  matrices 
are  functions  of  time,  or  time  varying. 

The  state  model  in  Equations  5.1  and  5.2  can  be  simplified  by  assuming  the 
plant,  control,  output,  and  feed-forward  matrices  are  not  functions  of  time.  By  using 
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this  assumption  the  Linear  Time  Invariant  (LTV)  state  model  can  be  described  as 


x(t)  =  Ax(t)  +  Bu(t ) 

(5.3) 

y(t)  =  Cx{t)  +  Du(t )  . 

(5.4) 

As  a  note,  the  state  model  representation  seen  in  Equations  5.3  and  5.4  is  the  most 
prevalent  matrix  form  encountered  when  dealing  with  linear  models,  as  most  systems 
are  simplified  by  assuming  time  invariance.  This  is  seen  in  modern  linear  optimal 
control  theory  [25,26]  and  optimal  state  estimation  theory  [9]. 

In  order  to  identify  the  state  vector  x(t)  Vi  >  0  for  linear  time  invariant  systems 
the  input  vector  u(t)  and  the  initial  state  vector  at  time  zero,  x(0 )  is  used,  as  seen 
in  Equation  5.5 

x(t)  =  $(i,  i0)x(0)  +  f  $(i,  t)Bu(t)cIt  .  (5.5) 

Jo 

In  the  above  Equation  5.5  the  State  Transition  matrix ($(i,i0))  is  required.  The 
state  transition  matrix  is  defined  for  a  linear  time  invariant  system  as  the  matrix 
exponential  in  Equation  5.6 


$(t)  =  eAt  (5.6) 

where  the  matrix  exponential  can  be  described  by  the  power  series  of  Equation  5.7. 

eAt  =  I  +  At  +  A2t 2  +  ^AV  +  •••  .  (5.7) 

Thus,  substituting  the  relation  depicted  in  Equation  5.6  into  Equation  5.5  the  expres¬ 
sion  for  the  state  vector  x(t)  is  then  redefined  as  seen  below. 

x(t)  =  eAtx(0)  +  f  eA(I~V Bu{t)(1t  (5.8) 

Jo 


44 


Using  the  same  procedure  as  the  state  vector,  the  system  output,  y(t),  is  defined  in 
a  similar  fashion. 


y(t)  =  CeAtx( 0)  +  [  CeA^  t-Bu(t)cIt  +  Du(t ) 


(5.9) 


Jo 


5.1.1  Transfer  Functions  for  LT1  Systems.  Determining  the  system  response 
for  LT1  systems  can  be  accomplished  by  way  of  the  method  outlined  in  Equation  5.9. 
This  method,  however,  operates  in  the  time  domain  and  requires  convolution,  which 
generally  limits  its  applicability.  By  transforming  the  LT1  system  into  the  frequency 
domain  by  the  Laplace  transform  (£)  the  process  of  determining  system  response  is 
greatly  simplified. 

In  order  to  determine  the  frequency  domain  transfer  function  of  the  system, 
G(s)  ,  the  Laplace  transform  of  the  time  domain  state  equations  in  Equations  5.3 
and  5.4  is  performed,  producing 


sX{s)  -  x(0)  =  AX(s)  +BU{s ) 
F(s)  =  CX(s)  +  DU(s)  . 


(5.10) 

(5.11) 


The  notation  X ,  Y  and  U  represent  the  Laplace  transforms  of  x,y  and  u.  With  the 
use  of  Equations  5.10  and  5.11  the  Laplace  transformed  output  Equation  5.9  is  then 
Y(s)  =  C[y(t)],  or  as  represented  below. 


y(s)  =  C{sl  -  A)-1x(0)  +  [C(sl  -  A)~1B  +  D]U(s)  (5.12) 


The  Laplace  transform  of  the  time  domain  impulse  response  matrix(g,(t))  results  in  the 
Laplace  transfer  function,  or  more  commonly  referred  to,  the  Transfer  function  (G( s)) 


G(s)  =  C[g(t)\ 

=  C(sl  -  A)~lB  +  D  . 


(5.13) 

(5.14) 
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By  using  the  relation  in  Equation  5.14  the  output  equation  Y(s)  in  Equation  5.12  can 
be  restated  as  seen  below. 


y(s)  =  G(s)U(s)  .  (5.15) 

5.1.2  Frequency  Response  for  LTI  Systems.  The  term  frequency  response 
refers  to  the  steady  state  response  of  a  system  to  a  sinusoidal  input  signal.  Frequency 
response  methods,  such  as  the  Bode  [25]  plot  display  the  system  response  over  a 
varying  range  of  input  frequencies.  The  system  response  is  generally  measured  in 
terms  of  the  magnitude  of  the  output  signal  and  the  phase  difference  between  the 
input  and  output. 

Frequency  response  characteristics  have  varying  properties  depending  on  the 
description  of  the  system  model.  In  the  case  of  LTI  systems,  the  input  and  output 
frequencies  are  the  same  while  having  different  magnitude  and  phase  differences  as 
noted  above.  Other  systems,  such  as  LTP  models,  the  input  and  output  frequencies 
vary.  This  is  an  important  distinction  and  will  be  covered  in 

For  LTI  systems  the  input  signal  used  for  frequency  response  is  the  exponentially 
modulated  sinus  oid{u0est)  as  seen  in  Equation  5.16. 

u(t)  =  u0est  s  G  C,  u0  E  Cm  .  (5.16) 

By  using  the  definition  of  the  Laplace  transfer  function  in  Equation  5.14  the  system 
output  can  be  described  as  the  product  of  it  and  the  input  of  Equation  5.16. 

Y(s)  =  G(s)u0est  .  (5-17) 

The  magnitude  for  single  input- single  output ( SISO)  systems  is  defined  as  the  modulus 
of  the  transfer  function  |G(s)|. 

|G(s)|  =  a Jreal(G(s))2  +  imag(G(s ))2  (5.18) 
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The  magnitude  of  the  system  is  also  referred  to  as  the  system  gain.  The  phase 
angle (Z(G(s)))  for  SISO  systems  is  defined  as  seen  in  Equation  5.19. 


Z(GW)  =  Arctan(‘^§m  (5.19) 

\  real(G(s))  J 

For  multi  input-multi  output  (MIMO)  systems  the  concept  of  magnitude  and  phase 
are  interpreted  differently  then  those  of  SISO  systems.  The  response  magnitude  for 
MIMO  system  is  defined  as  seen  in  5.20. 

Magnitude  =  -  f  (5.20) 

Nil 

where  ||  ■  ||  is  the  2-Norm  or  Euclidean  norm 

1 1 ex 1 1  =  yj~a\  +  .  (5.21) 

The  MIMO  system  response  is  different  from  that  of  the  SISO  system  as  the 
magnitude  response  is  defined  over  a  range 

minu o  <  Magnitude  <  maxUo  (5.22) 

as  the  magnitude  in  this  case  is  dependent  on  both  u0  and  {s  :  s  =  jco}.  Conven¬ 
tionally,  Uq  is  normalized  to  unity  ||m0||  =  1.  An  example  of  the  principal  gains  of  a 
MIMO  system  is  seen  seen  in  Figure  5.1  The  SVD  approach  outlined  in  Section  4.3  is 
commonly  used  as  a  method  for  determining  the  MIMO  system  gains  as  the  method 
outlined  in  Equation  5.20  can  be  overly  complex  for  a  large  systems.  By  performing 
the  SVD  on  the  transfer  function  G(s)  for  ever  frequency  over  {s  :  s  =  jut},  the  system 
principal  gains  are  then  the  frequency  dependent  singular  values,  a.  The  principal 
gains  are  determined  over  a  range,  as  outlined  in  Equation  5.22.  The  MIMO  phase 
relation  for  frequency  response  is  undefined  as  there  is  no  convenient  way  to  describe 
it.  This  is  due  to  each  input-output  relationship  having  a  different  phase  shift. 
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Singular  Values 


Figure  5.1:  Principal  Gains  of  a  MIMO  LTI  system  [2], 

5.2  Signal  Spaces  for  LTI  and  LTP  Systems 

Before  the  transfer  function  and  frequency  response  of  the  LTP  system  are 
explained  it  is  important  to  first  describe  the  differences  in  the  signal  spaces  of  LTI 
vs.  LTP  systems. 

A  LTI  system  having  a  sinusoidal  input  signal  of  frequency  (uf)  will  map  an 
output  signal  of  cun  but  with  different  magnitude  and  phase  then  the  input  signal. 
Thus,  the  signal  spaces  are  identical,  as  the  input  and  output  are  both  complex 
exponential  sinusoids  of  frequency  a if.  This  response  for  a  LTI  system  is  clearly  seen 
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in  a  Bode  plot,  as  phase  and  magnitude  of  the  system  are  plotted  over  the  entire 
range  of  input  frequency  ,  Uf  G  (0,  oo). 

A  LTP  system  differs  from  a  LT1  system  by  the  system  parameters  varying  pe¬ 
riodically  rather  than  remaining  constant,  or  time  invariant.  Consider  a  LTP  system 
defined  as  such  by  having  a  time  periodic  input  matrix,  B.  In  this  case  the  steady  state 
forced  response  signal  of  a  LTP  system  will  be  composed  of  the  sinusoidal  input  signal 
tOf  and  the  pumping  frequency  input  signal  c vp,  as  seen  in  Figure  5.2.  In  this  Figure 


1— 2cos/?o )  f 


Figure  5.2:  A  Simple  LTP  System  [42], 

an  input  signal,  u(t),  is  convoluted  with  a  time  periodic  gain  1  —  2j3cosu>pt.  Thus, 
the  output  signal,  y(t),  is  a  convolution  of  signals  best  described  by  Fourier  series 
containing  the  frequencies  Uf  ±  ncup  :  ujf  €  C,  n  G  Z,  which  is  an  infinite  dimension 
signal  space.  Clearly  the  input  and  output  spaces  are  not  the  same  dimension  as  in 
the  case  of  the  LTI  system  spaces.  This  difference  in  signal  spaces  differentiates  the 
LTP  system  from  LTI  systems,  as  the  mapping  of  signal  spaces  from  input  to  output 
is  not  an  isomorphism ,  being  one-to-one  and  onto,  but  rather  a  one-to-many  mapping 
as  seen  in  Figure  3.1.  To  overcome  the  difference  in  signal  spaces,  Wereley  [42]  devel- 
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oped  the  concepts  of  Geometrically  periodic  signals  and  the  identical  Exponentially 
modulated  periodic  signal ,  which  will  now  be  covered. 


5.2.1  Geometrically  Periodic  Signals.  From  Floqnet  analysis  of  a  LTP 
system  the  steady  state  response  at  the  time  t  is  related  to  the  system  response  a  full 
period  later  at  time  t  +  T  by  the  relation  in  Equation  5.23 

x(t  +  T)  =  zx(t),zE  C  (5.23) 

where  the  complex  scalar  z  G  C  is  defined  as 

2  =  esT  (5.24) 

and  as  a  note  s  =  juj.  By  using  the  relation  defined  in  Equation  5.23  the  system 
response  N  periods  later,  (t  +  NT),  is  defined  in  Equation  5.25 

x(t  +  T)  =  zNx(t),zeC.  (5.25) 

In  order  to  achieve  a  linear  mapping  Wereley  [42]  noted  the  input  signal  must  be  of 
the  form  of  Equation  5.25,  which  led  to  his  development  of  the  Geometrically  periodic 
signal,  as  defined  in  5.2.1 

Definition  5.2.1  (Geometrically  Periodic  Signals  [42]).  A  geometrically  periodic  sig¬ 
nal  (GP),  u(t),  having  a  pumping  frequency  ujp  and  period  T  is  defined  as 

u(t  +  T)  =  zNu{t),z  G  C  (5.26) 

if  there  exists  a  nonzero  z  G  C  and  N  e  N. 

5.2.2  Exponentially  Modulated  Periodic  Signals.  Wereley  [42]  further  pro¬ 
vided  that  the  definition  of  the  geometrically  periodic  signal  can  be  similarly  expressed 
as  a  complex  exponential  modulation  of  a  periodic  signal,  as  seen  in  definition  5.2.2 
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Definition  5.2.2  (Exponentially  Modulated  Periodic  Signals  [42]).  A  exponentially 
periodic  signal  (EMP),  u(t),  is  defined  as  the  complex  Fourier  series  expansion  of  a 
periodic  signal  having  a  pumping  frequency  oop  and  period  T,  modulated  by  a  complex 
exponential  signal  esTV  s  G  C 

OO 

u(t)  =  esT  unejnUpt 

n=— oo 

OO 

=  uneSnt  \t  >  0,  sn  =  s  +  jnujp  .  (5-27) 

n=—o o 

The  exponentially  modulated  periodic  signal  will  be  used  extensively  in  the 
development  of  the  the  state  space  form  of  the  LTP  system,  which  will  be  covered 
next. 

5.3  LTP  Transfer  Functions 

The  linear  operator  for  LTI  systems  is  the  well  known  Transfer  function  ,  G(s). 
As  was  covered  in  Section  5.1.1  the  Transfer  function  can  be  described  by  relating  the 
matrices  A,B,C,D  of  the  state  space  form  as  seen  in  Equation  5.14,  which  is  restated 
below 


G(s)  =  C(sl  -  A)_1B  +  D  . 

For  LTI  systems  the  frequency  response  is  directly  related  to  the  transfer  function 
G(s)  due  to  the  fact  that  in  steady  state  the  input  to  output  signal  maintains  the 
same  frequency  but  different  magnitude  and  phase.  Therefore,  the  input  and  output 
signal  spaces  are  of  the  same  dimension.  As  noted  in  Section  5.2  the  differences  in 
input  to  output  signal  spaces  precludes  the  direct  application  of  the  LTI  form  of  the 
transfer  function  for  LTP  systems.  By  modifying  the  input  signal  space  by  use  of 
the  EMP  signal  from  Section  5.2.2,  however,  a  convenient  form  similar  to  Equation 
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5.14  can  be  used  for  LTP  systems.  This  is  accomplished  by  the  method  of  Harmonic 
Balance 

5.3.1  Harmonic  Balance.  Analysis  techniques,  such  as  stability  analysis,  for 
LTP  systems  is  the  method  of  Harmonic  Balance.  This  method  originated  in  the  late 
19th  century  by  Hill  [11]  on  his  work  to  determine  the  stability  analysis  techniques  of 
linear  periodic  systems,  in  this  case  the  Lunar  perigee.  This  technique  will  be  briefly 
described  as  it  is  the  fundamental  building  block  in  the  development  of  the  state  space 
linear  operator  form  for  LTP  systems. 

Before  the  principle  of  Harmonic  Balance  is  defined  a  review  of  the  properties 
of  the  Fourier  series  must  first  be  covered.  Fourier  series  expansions,  being  either 
Trigonometric  or  complex,  of  harmonic  functions  involve  the  expansion  of  the  para¬ 
metric  excitation  into  a  set  of  basis  function.  These  basis  functions,  in  this  case 
sinusoids,  form  an  orthonormal  basis  over  the  fundamental  period.  As  was  covered 
in  Section  4.1.1,  the  elements  of  an  orthonormal  basis  are  linearly  independent  and 
therefore  the  coefficients  that  multiply  each  basis  function  in  a  Fourier  series  must 
vanish  for  each  harmonic  independently.  For  any  linear  periodic  system  whose  peri¬ 
odic  parameters  are  expanded  by  the  Fourier  series,  that  series  expansion  is  referred 
to  as  the  Harmonic  Balance.  The  Principle  of  Harmonic  Balance  refers  to  the  re¬ 
quirement  that  the  individual  coefficients  that  multiply  the  elements  of  the  Harmonic 
must  vanish  independently. 

The  principle  of  harmonic  balance  as  formed  by  Hill  [11]  used  a  exponential 
Fourier  series  which  was  modulated  by  a  complex  exponential.  This  can  be  seen  by 
determining  the  state  vector  solution  to  a  periodic  dynamic  equation,  such  as  the  Hill 
differential  equation  in  Equation  5.28 

x  T  [a  —  2 qip(t)}x(t)  =  0  (5.28) 
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where  ip{t)  is  defined  as  the  parametric  excitation{if{t ))  having  fundamental  period 
T.  Thus,  ip(t)  =  if(t  +  T).  Additionally,  the  parameter  q  is  referred  to  as  the 
pumping  amplitude  (q)  and  the  parameter  a  is  referred  to  as  a  constant  portion  of  the 
time  periodic  coefficient  of  x(t).  As  if{t)  is  T-Periodic,  it  can  be  described  as  by  the 
complex  Fourier  series,  as  seen  in  Equation  5.29. 

OO 

if{t)  =  J2  nejn“pt  (5.29) 

71= —  OO 

where  u ;p  is  the  pumping  frequency  of  the  harmonic  function.  The  assumed  solution  to 
the  periodic  Equation  5.28,  as  proposed  by  Hill  [11]  uses  the  Fourier  series  in  Equation 
5.30. 

OO 

x{t)  =  est  xnejnuJpt  ■  (5.30) 

71=  —  OO 

By  substituting  Equation  5.30  and  5.29  into  Equation  5.28  ,  and  noting  that  the  set 
of  functions  {ejnulp  :  n  G  Z}  forms  a  set  of  orthonormal  basis  functions  in  Lg  on 
the  interval  [ 0,7 ],  the  principal  of  harmonic  balance  can  be  applied  to  the  resulting 
equation,  as  seen  in  Equation  5.31 

OO 

0  =  [(s  +  jnup )2  +  a]xn  -  2q  ^  ^ mXn-m  ■  (5-31) 

771=  —  OO 

By  the  application  of  the  principle  of  harmonic  balance,  the  above  equation  results  in 
an  infinite  set  of  simultaneous  equations.  Hill  developed  the  determinant  of  the  infinite 
set  in  Equation  5.31  to  determine  the  stability  characteristics  of  periodic  functions 
as  a  function  of  the  periodic  parameters  of  the  function,  as  in  a  and  q  in  the  Hill 
equation. 

The  form  of  the  series  expansion  used  in  Equation  5.31  was  adopted  by  Wereley 
[42]  to  develop  a  linear  operator  for  periodic  systems  in  state  space  form.  This  was 
principally  done  by  adopting  the  use  of  the  exponentially  modulated  periodic  signal 
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from  Section  5.2.2  in  a  similar  manner  as  Hill  took  in  Equation  5.30.  This  will  be 
covered  next  as  the  Harmonic  Balance  State  Space  Model  for  LTP  systems  is  discussed. 

5.3.2  Harmonic  Balance  State  Space  Model.  As  discussed  in  Section  5.2 
the  signal  spaces  of  LTP  systems  are  not  equal,  which  precludes  the  development  of 
a  linear  operator  describing  the  input-output  mapping  of  a  signal  having  the  same 
input  and  output  spaces.  Wereley  [42]  has  demonstrated  that  by  using  a  exponentially 
modulated  periodic  signal  (EMP)  as  an  input  signal  to  a  LTP  system,  an  input  signal 
to  a  LTP  system  maps  a  EMP  input  signal  to  an  EMP  output  signal  having  the 
same  frequency  but  not  the  same  magnitude  or  phase  as  the  input  signal.  This 
is  significant,  as  the  input  and  output  signal  spaces  are  equal,  thus  the  concept  of 
a  linear  operator  is  now  possible  to  describe  for  a  LTP  system.  As  a  reminder, 
the  exponentially  modulated  periodic  signal  is  described  in  Section  5.2.2.  Also,  the 
magnitude  differences  in  input-output  mapping  of  a  LTP  system  refer  to  the  amplitude 
of  all  included  harmonics  in  the  input  and  output  signals.  As  a  linear  operator  has 
been  shown  to  exist  for  LTP  systems,  as  stated  above,  the  principal  of  harmonic 
balance  will  now  be  applied  to  transform  a  LTP  system  to  state  space  form. 

The  principal  element  to  provide  equal  input  and  output  signal  spaces  is  the 
EMP  input  signal,  u(t),  of  Equation  5.27,  and  restated  below 

OO 

u(t)  =  esT  unejnu)^ 

n=— oo 

OO 

=  uneSnt\  t  >  0, sn  =  s  +  jniVp  . 

TL— — OO 

Next,  by  using  an  approach  similar  to  Hill  to  provide  steady  state  response  to  of  a 
LTP  signal  to  an  EMP  signal,  as  in  Equation  5.30,  Wereley  [42]  and  later  Hwang  [12] 
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provide  the  steady  state  response  as 


x (t)  =  ^2  xneSnt  \t  >  0,  sn  =  s  +  jnujp  (5.32) 

n=— oo 
oo 

x(t)  =  ^2  snxneSrit  1 1  >  0 ,sn  —  s  +  jnujp  .  (5.33) 

n=— oo 

In  a  likewise  approach,  the  output  signal  of  a  LTP  system,  y(t),  is  a  linear  combination 
of  the  state  and  control  and  thus  is  described  as 

OO 

y(t)  =  ^2  yneSnt  \t>0,sn  =  s  +  jnujp  .  (5.34) 

n=— oo 

Now,  considering  the  form  of  a  LTP  system  in  state  space  form  in  Equations  5.1  and 
5.2,  as  restated  below 

x(t)  =  A(t)x{t)  +  B(t)u{t) 
y{t)  =  C(t)x(t)  +  D(t)u(t )  . 

It  is  important  to  note  that  the  matrices  [A(t),  B(t),  C(t),  and  D(t)]  are  T-periodic 
and  can  be  expanded  by  the  Fourier  series,  as  seen  in  Equation  5.35  for  the  plant 
matrix,  A(t). 

OO 

A(t)  =  "^2  \t  >  0,  sn  =  s  +  jnujp  .  (5.35) 

n=— oo 

As  a  note,  all  the  remaining  matrices  [B(t),  C(t),  and  D(t)]  can  be  expressed  as  in 
Equation  5.35.  Next,  by  substituting  Equations  5.32,  5.33,  5.34  and  5.35  along  with 
the  similar  approaches  for  [B(t),  C(t),  and  D(t)]  the  LTP  state  space  form  can  be 
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restated  as  in  Equation  5.36 


E 


rp  psnt  - 


Smt 


Anesnt  XmeSmt  +  Y1  BneSnt  Y1 

n=— oo  m=— oo  n=— oo  m=— oo 

oo  oo 

E/1  zy>  I  \  n ,  pSn-|-mt 

-rx72*-t/772c'  i  /  J-JnUjm'z' 

n=— oo  72,722=— oo 

oo  oo 

An_mxmeSnt  +  Bn_mumeSnt  .  (5.36) 


72,772=— oo 
oo 


n,m=— oo 


n,m=—o o 


A  likewise  approach  can  be  made  for  the  output,  y(t),  as  seen  in  Equation  5.37 


n  —  En,m=-o o  Cn-mXm,eSnt  +  En,m=- oo  Dn_mUmeSnt  .  (5.37) 


Next,  by  moving  all  terms  of  Equation  5.36  and  5.37  to  the  right  hand  side  and  then 
multiplying  through  by  e_st  results  in 


0 

0 


E 

7.—  —  OC 
OO 

E 


OO  oo 

Sn%n  ^  ^  A n—m%m  “t-  ^  ^ 

m=— oo  m=— oo 

oo  oo 

Vn  ^  ^  C'n—mXm  +  ^  ^  D n_rriUrn 

m=—oo  m=— oo 


Snt 


(5.38) 

(5.39) 


Now,  by  the  principle  of  harmonic  balance  from  Section  5.31,  the  terms  within  the 
brackets  in  Equation  5.39  are  linearly  independent  and  therefore  must  be  null  to 
avoid  a  trivial  solution.  Therefore,  the  state  model  of  Equation  5.39  can  be  restated 
Vn  G  Z  as 


OO  OO 


^  ^  -4-n— + 

772—  —  OO 

^  ^  4^72—772^772 

772=  —  OO 

(5.40) 

OO 

^  ^  ^72  —  772*^772  4" 

OO 

^  ^  -^72—772^772  * 

(5.41) 

m=— oo  m=— oo 


The  input-output  relation  stated  in  Equation  5.41,  while  functional,  is  difficult  to  use 
due  to  the  summations.  By  applying  the  Toeplitz  transformation  from  Section  4.4  to 
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the  T-periodic  matrices  A(t),  B(t),  C(t),  D(t ),  Wereley  [42]  simplified  Equation  5.41, 
as  defined  below. 

Definition  5.3.1  (Harmonic  State  Space  Model  [42]).  The  Toeplitz  transformation  of 
the  T-Periodic  matrices  of  Equation  5.fl  are  expressed  as  the  doubly  infinite  matrix 
equation 


sx  =  (A  —  Af)x  +  Bu  (5.42) 

y  =  Cx  +  Vu  (5.43) 

where  J\f ,  originally  defined  in  Equation  f.26,  represents  a  doubly  infinite  block  diag¬ 
onal  matrix  consisting  of  all  of  the  pumping  frequency  harmonics,  ncop 

AT  =  blkdiag{jnujpI}\/ n  <E  Z  .  (5.44) 

Additionally  the  state,  control,  and  output  vectors  are  doubly  infinite,  as  stated  as 


X-2 

u- 2 

y  2 

X-i 

W-l 

y  i 

x0 

U  = 

u0 

y  = 

y0 

X1 

U\ 

2/i 

X2 

U-2 

2/2 

The  Toeplitz  transformed  T-periodic  A,B,C,D  matrices  are  made  of  the  complex  Fourier 
coefficients  of  the  matrices,  as  covered  in  Section  4.22.  For  example,  the  plant  matrix 
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Aft)  is  defined  as 


Ao 

A_! 

A_2 

1 

1 

CO 

1 

A ! 

Ao 

A_! 

1 

to 

1 

CO 

A'2 

Ax 

Ao 

A_x  A_2  •  •  • 

(5.45) 

A3 

A2 

Ax 

Ao 

A„x 

a4 

A3 

A2 

Ax 

A0 

where  the  matrices  B,C,and  D  are  described  using  the  same  methodology  as  described 
in  Equation  5-45. 

Now  that  the  harmonic  state  space  form  has  been  defined  in  a  manner  similar 
to  the  LT1  form,  the  transfer  function  for  periodic  systems  can  now  be  defined  in  a 
a  manner  similar  to  a  LTI  system.  By  using  the  LTP  harmonic  state  space  model 
from  definition  5.3.1  the  input-output  relationship  between  the  input  harmonics,  {un  : 
n  G  Z},  and  the  output  harmonics,  {yn  :  n  G  Z},  can  be  described  by  the  harmonic 
transfer  function,  G(s). 

Definition  5.3.2  (Harmonic  Transfer  Function).  The  harmonic  transfer  function  , 
G(s),  is  an  infinite  dimensional  matrix  of  Fourier  coefficients  relating  the  the  input 
harmonics,  {un  :  n  6  Z},  and  the  output  harmonics,  {yn  :  n  6  Z},  described  as 

y  =  G(s)u  (5.46) 

where 

G(s)  =  C[sl  -  (A  -  M)}~lB  +  V  (5.47) 

as  long  as  the  inverse  [si  —  (A  —  A/")]-1  exists. 
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This  chapter  presented  the  fundamentals  of  linear  time  periodic  systems,  with 
an  emphasis  on  the  state  space  form.  The  next  chapter  will  develop  the  Cramer- 
Rao  bound  in  a  manner  similar  to  LTI  systems  by  using  the  LTP  state  space  form 
developed  in  this  chapter. 
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VI.  Linear  Time  Periodic  System  Parameter  Validation  via 
the  Cramer-Rao  Lower  Bound 


The  intent  of  this  researdi,  as  stated  previously,  is  to  develop  a  more  accurate 
method  of  rotor  smoothing  by  way  of  controlling  a  verifiably  accurate  rotor 
system  model.  In  the  last  chapter  the  concept  of  the  linear  time  periodic  model  was 
developed,  which  provides  the  foundation  for  an  accurate  helicopter  rotor  in  forward 
flight.  The  identified  parameters  that  populate  that  LTP  model,  however,  must  be 
accurate  in  order  to  provide  a  basis  for  an  effective  controller.  This  chapter  will 
present  the  development  of  the  Cramer-Rao  lower  bound  for  a  linear  periodic  system 
in  order  to  validate  the  identified  system  parameters. 

6. 1  Introduction 

Helicopter  rotors  are  subject  to  vibratory  forces  from  a  variety  of  sources.  In 
forward  flight,  vibratory  forces  at  integer  harmonics  of  the  rotor  speed  are  generated 
because  the  velocity  of  the  blades  with  respect  to  the  air  is  inherently  periodic.  In 
hover,  vibratory  forces  in  the  plane  of  the  rotor  disk  are  generated  as  a  result  of  mass 
imbalances  in  the  rotor  blades;  and  out-of-plane  vibratory  forces  result  from  an  mi- 
symmetric  distribution  of  lift  caused  by  aerodynamic  dissimilarities  among  the  blades. 
These  vibratory  forces  have  a  period  equal  to  the  rotor  period  of  revolution  (one-per- 
rev).  While  other  sources  of  vibration  also  exist,  the  main  rotor  system  generates  by 
far  vibrations  of  the  largest  magnitude.  Over  time,  the  vibrations  produced  by  the 
main  rotor  will  result  in  damage  to  the  aircraft,  and  can  compromise  the  effectiveness 
of  the  crew.  While  vibratory  loads  in  forward  flight  are  to  some  degree  unavoid¬ 
able,  the  one-per-rev  vibratory  loads  due  to  mass  and  aerodynamic  imbalances  can 
be  minimized  through  proper  and  regular  maintenance. 

Improving  maintenance  procedures  for  eliminating  one-per-rev  vibrations  for 
helicopter  main  rotors  (rotor  smoothing)  has  been  an  area  of  active  interest  since 
the  early  days  of  helicopters,  when  those  procedures  were  called  track  and  balance. 
The  current  state-of-the-art  in  rotor  smoothing  techniques,  while  improving,  is  still 
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costly  and  time  consuming.  One  of  the  principle  reasons  for  the  inefficiency  of  current 
rotor  smoothing  maintenance  methods  is  that  they  all  require  multiple  iterations;  each 
necessitating  a  separate  test  flight  followed  by  maintenance  actions  to  adjust  blade 
mass,  pitch  link  length,  or  trim  tab  angle,  until  the  vibrations  are  reduced  to  an 
acceptable  level.  In  order  to  minimize  the  number  of  iterations  required  to  achieve  an 
acceptable  vibration  level,  one  must  accurately  predict  the  optimal  set  of  adjustments. 
These  predictions  rely  exclusively  on  a  having  an  accurate  plant  model  of  the  entire 
aircraft,  since  fuselage  characteristics  affect  where  and  how  much  vibration  is  sensed. 

In  order  to  develop  an  accurate  model  of  a  helicopter,  the  parameters  that  define 
the  system  must  be  accurately  identified.  While  it  is  obvious  that  different  helicopter 
types  require  different  models,  it  may  be  less  obvious  that  individual  helicopters 
of  the  same  type  may  have  significant  differences  that  require  modifications  to  the 
parameters,  and  that  those  parameters  may  change  over  time  [18].  Therefore,  it  is 
important  to  be  able  to  identify  an  aircraft’s  parameters  ‘on  the  fly’.  Several  of  the 
most  current  rotor  smoothing  methods  use  system  identification  methods  based  on 
either  a  linear  least  squares  approximation  [24],  or  a  neural  network  approximation  [39, 
43].  Due  to  the  high  cost  of  flight  tests,  these  methods  typically  depend  on  a  very 
small  dataset,  usually  obtained  from  a  single  aircraft. 

Regardless  of  the  size  of  the  dataset  used  to  create  the  aircraft  model,  or  the 
number  of  different  aircraft  used,  it  is  equally  important  that  the  parameters  identified 
for  an  aircraft  are  verihably  accurate.  For  linear,  time  invariant  (LTI)  systems,  the 
accuracy  of  identified  parameters  can  be  verified  using  the  Cramer-Rao  bound.  How¬ 
ever,  a  helicopter  rotor  is  time  periodic,  so  the  Cramer-Rao  bound  for  LTI  systems  is 
not  valid.  The  objective  of  this  chapter  is  to  develop  the  theory  necessary  to  define 
the  Cramer-Rao  bound  for  a  linear,  time  periodic  (LTP)  system.  With  such  a  tool, 
plant  models  derived  from  identified  parameters  for  LTP  systems  can  be  verified  by 
examining  the  Cramer-Rao  bound  of  each  parameter  for  accuracy,  and  adjustments 
can  be  made  to  improve  model  accuracy. 
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This  chapter  begins  briefly  describing  the  Cramer-Rao  bound  as  used  to  deter¬ 
mine  parameter  accuracy.  Next,  the  Maximum,  Likelihood  Estimator  (MILE)  will  be 
defined,  as  it  forms  a  direct  link  to  the  formulation  of  the  Cramer-Rao  lower  bound. 
Next,  to  provide  a  solid  foundation  for  the  reader,  the  definitions  of  several  statis¬ 
tical  methods  used  to  define  the  Cramer-Rao  inequality  are  presented.  In  turn,  the 
Cramer-Rao  inequality  is  then  formulated  based  on  the  MLE.  Finally,  by  way  of  the 
definition  of  the  Cramer-Rao  inequality,  the  Cramer-Rao  lower  bound  is  defined  for 
individual  parameters  of  a  LTP  system  in  state  space  form.  The  effectiveness  of  the 
Cramer-Rao  lower  bound  is  demonstrated  at  the  conclusion  to  this  chapter  through 
a  short  example  depicting  parameter  accuracy  for  a  LTP  system. 

6.2  General  Description  of  the  Cramer-Rao  Lower  Bound 

The  parameters  determined  by  system  identification  methods  define  a  linearized 
vehicle  model,  which  in  turn  defines  the  dynamic  characteristics  of  that  vehicle.  These 
models  are  then  used  in  simulators,  control  system  design,  and  are  used  to  validate 
wind  tunnel  parameter  predictions.  Maine  and  Iliff  [34]  point  out  that  it  is  impor¬ 
tant  to  remember  that  parameters  obtained  from  testing  are  only  estimates,  and  not 
exact  values.  This  fact,  unfortunately,  is  often  disregarded  and  only  rarely  are  the 
parameters  ever  verified  for  accuracy.  Maine  and  Iliff  [34]  further  state  that  if  accu¬ 
rate  parameter  estimates  cannot  be  distinguished  from  worthless  estimates,  then  all 
estimates  must  be  assumed  to  be  of  questionable  accuracy.  It  is  for  these  reasons  that 
parameter  validation  methods  have  been  developed  for  linear,  time  invariant  (LTI) 
systems  such  as  those  cited  in  Refs.  [34]  and  [3]. 

As  stated  in  Section  6.1,  the  intent  of  this  research  is  to  improve  LTP  modeling 
for  controller  development,  with  specific  application  for  alleviating  main  rotor  vibra¬ 
tions  through  the  use  of  rotor  smoothing  techniques.  The  measure  of  accuracy  of 
the  model  parameters  of  the  identified  rotor  system  are  of  critical  importance  in  this 
application  as  they  will  determine  the  effectiveness  of  the  vibration  control  implemen¬ 
tation.  Implementing  a  controller  based  on  an  unverified  model  may  have  disastrous 
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results,  as  the  commanded  control  inputs  are  based  on  a  model  that  does  not  match 
the  true  rotor  system  dynamics. 

The  Cramer-Rao  bound  is  a  method  commonly  used  in  flight  testing  to  establish 
the  accuracy  of  identified  parameters  of  a  linearized  vehicle  model.  As  demonstrated 
by  Maine  and  Iliff  [34],  the  Cramer-Rao  lower  bound  gives  the  theoretical  lower  limit 
to  the  accuracy  of  parameter  estimates  by  an  optimal  estimator.  This  measure  of 
accuracy  is  based  on  the  uncertainty  ellipsoid  and  is  similar  to  other  measures  of 
accuracy  such  as  estimated  variance  and  standard  deviation.  The  Cramer-Rao  bound 
is  similar  to  these  methods  except  that  the  Cramer-Rao  bound  is  the  square  root  of  the 
variance.  By  examining  the  Cramer-Rao  bound  for  each  parameter  estimated  from 
flight  test,  one  can  determine  the  accuracy  of  the  data  by  evaluating  the  size  of  the 
bounds.  The  Cramer-Rao  bound  establishes  the  standard  deviation  of  an  identified 
model  parameter,  therefore,  large  bounds  indicate  poor  estimation  performance.  This 
method  provides  greater  insight  into  the  accuracy  of  estimated  parameters  than  simply 
examining  the  scatter  of  parameter  estimates  from  flight  test  data  at  each  point,  as 
seen  in  Figure  6.1. 

Upon  review  of  Figure  6.1,  one  can  see  that  the  Cramer-Rao  bound  indicates 
that  data  below  -10  degrees  is  unreliable,  even  though  the  data  scatter  is  not  very 
large.  Thus,  an  examination  of  the  Cramer-Rao  bounds  for  estimated  vehicle  model 
parameters  provides  a  basis  for  an  evaluation  of  whether  the  parameters  have  been 
estimated  properly,  or  whether  they  should  either  be  re-evaluated  or  adjusted  by  the 
engineer  before  they  are  used  in  the  vibration  control  model. 

6.3  The  Maximum  Likelihood  Estimator 

As  stated  by  References  [34]  and  [10],  the  Maximum  Likelihood  Estimator (MLE) 
is  by  far  the  most  commonly  used  technique  for  estimating  parameters.  The  MLE  is 
a  method  that  relates  the  a  sampled  vector  of  Independent  arid  Identically  Distributed 
random  values(iid)  to  a  vector  of  ’true’  values  in  an  attempt  to  discern  differences 
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Solid  symbols  indicate  derivatives 
obtained  immediately  before 
aircraft  departures 


O  Flight  data 

Cram^r-Rao  bound 
-  Fairing 

Dashes  indicate  less  certain  fairings 


Figure  6.1:  Example  of  Cramer-Rao  Bounds  for  Parameter  Estimates  [34], 

between  the  two  vectors.  Before  the  MLE  is  defined,  the  concept  of  an  iid  sample  is 
first  defined  as  follows: 

Definition  6.3.1  (Independent  and  indentically  distributed  random  variables(iid)). 
The  random  variables  Xu ... ,  Xn  are  defined  as  independent  and  identically  dis¬ 
tributed  random  variables  if  Xi, ...  ,Xn  are  mutually  independent  random  variables 
from  a  probability  density  function(pdf)  of  the  same  function  F(x). 

Thus,  with  a  iid  sample  Xi, . . . ,  Xn  from  a  pdf  f(x\6i, . . . ,  9k),  the  likelihood 
function  is  defined  according  to  reference  [10]  in  Equation  6.1. 

n 

l(9 |x)  =  Y[f(xi\e1,...,ek) .  (6.i) 

2=1 

The  likelihood  function  defined  in  Equation  6.1  is  simply  a  measure  of  the  plausibility, 
or  likelihood ,  of  the  variable  0t  to  the  vector  of  sample  data  points,  X.  Having  defined 
the  likelihood  function,  the  maximum  likelihood  estimator  can  now  be  stated  as  in 
definition  6.3.2. 
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Definition  6.3.2  (Maximum  Likelihood  Estimator).  The  maximum  likelihood  esti¬ 
mator  of  the  variable  6  is  defined  as  the  estimate  of  the  value  9i, . . . ,  9k  that  maximises 
the  likelihood  function, L(9\x) .  Thus,  for  each  sample  point  x,  let  9(x)  be  a  parameter 
value  at  which  L{9 \x)  attains  a  maximum  value  as  a  function  of  9  while  holding  x 
fixed. 

The  MLE  can  be  produced  by  differentiation  of  the  likelihood  function  ,  however  direct 
maximization  through  an  optimization  function  is  generally  the  method  of  choice  to 
identify  the  parameter  6. 

The  form  of  the  likelihood  function  commonly  used  for  the  MLE  is  based  on 
Gaussian  distribution,  and  is  referred  to  as  the  normal  likelihood  function.  For  the 
construction  of  the  normal  likelihood  function,  let  x  —  9(x)  be  the  predicted  output 
based  on  the  postulated  value  6.  Thus,  the  normal  likelihood  function  is  defined  as 


L(d|x) 


Y[f(xi\9i,...,9k) 

2=1 


n 

2=1 


1  l/2 


2n\R\ 

n/2 


2tt\R\ 


e-|  [(£i-Xi)TR  Rxj— ®j)] 
e~\YG:=i[(£i-xi)T  R^Gi-xi)] 


(6.2) 


where  R  denotes  the  prediction  error  covariance  matrix ( R)  and  T  denotes  the  trans¬ 
pose  operator.  As  a  note,  the  prediction  error  covariance  matrix  plays  a  major  part 
in  the  Cramer-Rao  bound  and  is  covered  in  more  detail  in  Section  6.5.  Reference  [34] 
indicates  that  if  the  prediction  error  covariance  matrix  is  known,  the  term  multiply¬ 
ing  the  exponential  in  Equation  6.2  is  a  constant  and  therefore  can  be  discarded  as  it 
does  not  effect  the  maximization.  The  normal  likelihood  Equation,  6.2,  is  commonly 
modified  by  taking  the  negative  of  the  logarithm  to  create  the  logarithmic  likelihood 
function ,  as  seen  in  as  seen  in  Equation  6.3.  By  using  the  log  likelihood  function  as 
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a  cost  functional,  J(9),  the  MLE  can  be  obtained  through  direct  minimization. 

1  n 

m  =  £  [(f«  -  xi)TR~i{xi  -  a?,)]  (6.3) 

2—1 

The  cost  functional  stated  in  Equation  6.3  provides  a  direct  link  to  the  Cramer-Rao 
inequality,  as  will  be  seen  in  Equation  6.19.  While  the  MLE  provides  the  most  likely 
estimate  of  a  system  parameter,  it  provides  no  method  to  validate  the  accuracy  of 
that  estimate.  This  will  be  addressed  in  the  next  Section  regarding  the  Cramer-Rao 
inequality. 


6-4  The  Cramer-Rao  Inequality 

The  Cramer-Rao  inequality  is  a  statistical  method  which  provides  a  means  to 
evaluate  the  accuracy  of  a  given  parameter  estimate.  As  stated  by  reference  [34],  the 
Cramer-Rao  inequality  gives  a  theoretical  limit  for  the  accuracy  that  is  possible  for 
the  estimate  regardless  of  the  estimator  used.  Before  the  Cramer-Rao  inequality  is 
described,  the  concept  of  Expected  Values  ( Ex ,  hx),  Variance  (Var  A,  cr^-),  and  Co- 
variance(Cov(X,  YJ)  must  first  be  defined  as  they  are  used  extensively  in  this  Section. 
As  a  note,  the  definitions  listed  below  are  adapted  from  reference  [10]. 


6-4-1  Statistical  Methods  Defined.  The  expected  value  or  mean  of  a  random 
variable  is  simply  the  average  value  that  would  be  expected  from  a  random  sample  of 
a  probability  distribution.  The  expected  value  is  also  known  as  the  first  moment  of 
a  distribution,  and  is  defined  below. 


Definition  6.4.1  (Expected  Value).  The  expected  value,  E[x(t)\,  of  a  random  variable 
x(t)  of  the  probability  distribution  fx{x\t )  is  defined  as 


Jf^x(t)fx(x;t)dx  for  continuousdistributions 
Thx&x  X(t)fx(x]  t)  for  discrete  distributions 
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The  expected  value  is  generally  denoted  by  Ex  or  /4,,as  will  be  seen  throughout  this 
work. 

The  second  moment  of  a  distribution  is  known  as  the  variance ,  and  is  described 
as  in  definition  6.4.2 

Definition  6.4.2  (Variance).  The  variance  of  a  random  variable  x{t)  is  defined  as 

Var[x(t)\  =  E[x(t )  —  E(x(t))]2 
=  E[x(t)  -  yx}2  . 


The  variance  of  a  random  variable  x  is  also  denoted  as  cr2.  The  variance  defines  the 
spread  of  a  distribution  about  the  mean  of  the  random  variable,  jix .  The  spread  about 
the  mean  is  also  commonly  measured  by  the  square  root  of  the  variance,  known  as 
the  standard  deviation ,  ax. 

When  considering  two  random  variables,  x(t)  and  y(t ),  the  relationship  between 
the  two  can  be  evaluated  by  use  of  the  covariance ,  as  defined  below  in  definition  6.4.3 

Definition  6.4.3  (Covariance).  The  covariance  of  two  random  variables  x(t)  and  y(t ) 
is  defined  as 


Cov[x(t),y(t )] 


E[(x(t)  -  E(x(t)))(y(t)  -  E(y(t)))] 
E[(x(t)  -  nx){y{t)  -  Hy)\  . 


The  sign  of  the  covariance  relates  the  relationship  between  the  two  random  variables. 
Thus,  if  large  values  of  the  random  variable  x  are  observed  with  small  values  of  the 
random  variable  y,  the  sign  of  Cov[x(t),y(t)]  will  be  negative.  Likewise,  if  large  values 
of  the  random  variable  x  are  observed  with  large  values  of  the  random  variable  y,  the 
sign  of  Cov[x(t),y(t)}  will  be  positive. 
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6-4-2  The  Cramer-Rao  Inequality.  Using  the  statical  methods  described  in 
Section  6.4.1,  the  Cramer-Rao  lower  bound  will  now  be  defined  based  on  the  presen¬ 
tation  in  [10]. 

Definition  6.4.4  (The  Cramer-Rao  Inequality).  Let  X\, ... ,  Xn  be  a  sample  with 
the  probability  distribution,  f(x\6),  and  let  W{X)  =  W(X i, . . .  ,Xn)  be  any  estimator 
satisfying 

IesW(X)  =  I^L[W(x)f(x\e)}dx  (6.4) 

and 


VargW(X )  <  oo 


(6.5) 


Then 


VareW(X )  > 


/ 

m&og  Jvmyj ' 


(6.6) 


Proof:  The  proof  of  the  Cramer-Rao  inequality  is  provided  by  Casella  and  Berger  [10] 
and  is  stated  below.  The  proof  starts  by  using  the  Cauchy-Schwartz  Inequality,  which 
is  stated  in  Equation  6.7  for  any  two  random  variables  x(t)  =  X ,  y(t)  =  Y 

[Cao{X,Y)f  <  (Var  X)(VarY)  .  (6.7) 

Equation  6.7  can  be  rearranged  to  define  a  lower  bound  on  the  variance  of  the  random 
variable,  X ,  as  stated  in  Equation  6.8. 

[CovjX,  Y)}2 
( Var  Y) 
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( VarX )  > 


(6.8) 


Before  expanding  on  Equation  6.8  by  assigning  X  to  be  the  estimator  W(X)  and  Y 
to  be  the  quantity  j^log  f(X\9),  we  first  note  on  the  covariance  between  W(X)  and 
■§e log  f(X\9)  below: 


TeE°w{x) 


ms(x  m 


dx 


(6.9) 


Equation  6.9  is  modified  by  multiplying  through  by 


fW0) 

fix  |0)> 


which  results  in 


Toe°w^ 


Ee 


W{X) 


|CT) 

me) 


(6.10) 


Next,  Equation  6.10  is  transformed  by  the  property  of  logs  to  reveal  Equation  6.11 


fgEeWm  = 


Eg 


W(X)ff0log  f(X\e) 


(6.11) 


To  declare  Equation  6.11  as  a  covariance  between  W (A")  and  -§glog  f(X\9)  the  product 
of  the  expected  values  must  be  subtracted.  By  applying  W (x)  =  [1]  to  the  Equation 
6.11  the  expectation  of  j^log  f(X\9)  can  be  calculated  as 


Eel&log  f(X\0)]  =  jgEtl  1]  =  0.  (6.12) 

Thus,  the  expectations  above  are  null,  and  therefore  the  covariance, 

Covg(W(X),  ^ log  f(X\9)),  is  equal  to  the  expectation  seen  in  Equation  6.11.  This 
relationship  is  depicted  in  Equation  6.13. 


Cm,(W(X),  f  leg  f(X\»))  =  E,  [W(X)&log  f(X\0)]  =  &E,W(X)  (6.13) 


Noting  that 


d 

d9 


EeW(X )  =  I 


(6.14) 
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where  /  is  the  identity  matrix^ I).  Additionally,  as  Eg  [-§glog  f(X\ 0)]  =  0  the  variance 
of  fjylog  f(X\9 )  is  described  as  seen  in  Equation  6.15. 


Vare^log  f(X |0))  = 


Eg 


pogf(Xle))2 


(6.15) 


Finally,  substituting  Equations  6.13  and  6.15  into  the  covariance  and  variance  on  the 
right  hand  Equation  6.8,  the  Cramer- Rao  inequality  is  defined  in  Equation  6.16. 


Varg(W(X))  > 
> 


(JjEewyqf 

Es  [(& log  /( X\e)Y] 

I 

E,  [(| log  /(.V|9))2] 


(6.16) 


It  is  important  to  note  that  while  this  derivation  of  the  Cramer-Rao  inequality  is 
presented  for  continuous  time  random  variables  it  is  also  valid  for  discrete  time  random 
variables.  This  is  due  to  Equation  6.4,  in  which  the  integration  can  be  substituted 
for  summation  in  the  discrete  time  case. 

Two  important  facts  can  now  be  discussed  concerning  the  Cramer-Rao  inequal¬ 
ity.  The  first  is  that  Equation  6.16  presents  the  lower  bound  for  the  variance  of  any 
estimator.  Furthermore,  any  unbiased  estimator  that  attains  the  equality  of  Equa¬ 
tion  6.16  is  deemed  an  efficient  estimator.  In  the  case  of  this  work,  the  estimator  of 
interest  will  be  the  maximum  likelihood  estimator.  This  condition  is  met  based  on 
definition  6.4.5 

Definition  6.4.5  (Existence  of  an  Efficient  Estimator).  If  an  efficient  estimator  exists 
for  a  problem,  that  estimator  is  a  maximum  likelihood  estimator 

Proof:  See  reference  [34] 

The  second  fact  of  importance  from  Equation  6.16  is  the  existence  of  the  Fisher 
Information  Matrix ( M)  in  the  Cramer-Rao  inequality.  The  Fisher  information  ma- 
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tri x,M,  is  defined  as 


Me(X) 


(i llogf(X\e)f 


(6.17) 


Reference  [10]  notes  the  Fisher  information  matrix  is  referred  to  as  the  Fisher  in¬ 
formation  of  the  sample,  as  the  Fisher  information  value  defines  the  bound  of  the 
variance  of  the  estimator  of  6.  Thus,  as  the  information  number  gets  bigger,  more 
information  is  available  about  the  6.  As  a  result,  a  large  Fisher  number  results  in 
a  smaller  bound  on  the  variance  and  therefore  indicates  a  more  accurate  estimate 
of  6.  Now,  the  Cramer-Rao  inequality  can  be  restated  in  terms  of  the  the  Fisher 
Information  Matrix  as 


Vare(W(X ))  >  Mo(X)-1  .  (6.18) 

These  results  of  the  Cramer-Rao  inequality  indicate  that  the  lower  bound  of  the 
variance  is  a  close  approximation  of  the  variance  of  the  estimates  from  a  maximum 
likelihood  estimator.  This  presents  an  ideal  form  to  used  to  evaluate  the  performance 
of  parameter  estimates,  however,  the  form  of  the  Fisher  information  matrix  can  be 
cumbersome  to  compute.  Reference  [34]  presents  a  close  approximate  for  the  Fisher 
information  matrix  by  way  of  the  Hessian,  H,  of  the  cost  function,  J,  of  Equation  6.3, 
presented  as 


H  =  V2eJ  .  (6.19) 

This  approximation  is  used  in  the  computation  of  the  Cramer-Rao  lower  bound  in 
Section  6.5. 

6.5  Cramer-Rao  Lower  Bound  for  LTI  Systems 

The  following  is  an  adaptation  from  Ref.  3  that  describes  how  the  Cramer-Rao 
bound  is  calculated  for  a  LTI  system  is  state  space  form.  For  LTI  systems,  the  first 
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step  in  generating  the  Cramer- Rao  bound  is  to  generate  the  frequency  response  of  the 
linear  system  represented  by  Eq.  6.20,  which  is  in  state  space  form. 

x  =  Ax(t)  +  Bu{t) 

y  =  Cx(t )  +  Du(t )  (6.20) 

ffere,  the  frequency  response  of  the  linear  system  is  defined  as 

G(uuf)  =  C{iujfI  -  A)-lB  +  D  (6.21) 

where  it  is  assumed  that  the  matrix  D  is  zero,  since  no  feed  forward  signal  is  assumed 
in  this  model.  Next,  we  define  in  Eq.  6.22  the  error  between  the  estimated  model 
frequency  responses,  and  the  test  frequency  responses,  Gj^co^Test,  from 

the  real  system. 


^  f)  Gj!k{tUf)  Gjtk(lMf)test  (6.22) 

Here  the  indices  j  and  k  represent  the  number  of  frequency  response  outputs  and 
frequency  response  inputs,  respectively,  to  the  frequency  response,  Gj^{iuif).  The 
response  error,  £jtk  is  required  to  establish  a  cost  function  defining  the  fitness  of  the 
estimated  frequency  response  as  compared  to  the  actual  system  response,  J,  as  defined 
in  Equation  6.23.  Note  that  nQ  is  the  number  of  response  outputs,  nL  is  the  number 
of  control  inputs,  and  nu}f  represents  the  number  of  test  frequency  points. 

na  nL  n“/ 

j  =  E  E  E  [%*(^)]2  (6-23) 

j= 1  k= 1  1=1 

The  cost  function  of  Equation  6.23  is  based  upon  the  negative  of  the  logarithm  of 
the  maximum  likelihood  functional,  as  developed  in  Maine  and  Iliff  [34],  This  is  due 
to  the  fact  that  the  Cramer-Rao  bound  produces  the  theoretical  lower  limit  of  the 
accuracy  of  an  optimal  estimator.  In  the  case  of  the  Cramer-Rao  bound,  the  optimal 
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estimator  is  the  maximum  likelihood  estimator.  As  in  the  case  of  the  maximum 
likelihood  estimator,  the  cost  functional  of  Equation  6.23  assumes  a  Gaussian  white 
noise  measurement  disturbance  having  zero  mean.  The  intensity  of  the  measurement 
noise  vector  is  described  by  the  prediction  error  covariance  matrix,  R,  which  is  formed 
by  setting  R  equal  to  the  power  spectral  density  of  the  measurement  noise,  Sv.  The 
cost  functional  of  Equation  6.23  provides  input  for  the  prediction  error  covariance 
matrix  by  the  weighting  matrix  W.  The  matrix  W  represents  a  weighting  of  the 
inverse  of  input  noise  strength,  which  is  formed  by  W  =  R~1. 

The  frequency  response  fitness  function  J,  defined  in  Eq.  6.23  forms  the  core  of 
the  Cramer-Rao  bound  by  providing  the  necessary  elements  for  the  Hessian,  H.  The 
Hessian  of  J  with  respect  to  a  selected  model  parameter  p  is  defined  in  Eq.  6.24. 


H  = 


(6.24) 


Thus,  the  Cramer-Rao  bound,  CRj,  is  approximated  in  terms  of  the  Hessian  for  each 
identified  system  parameter  of  Eq.  6.21,  pj  as 

CRj^  yJ(S-')u  ■  (6.25) 


In  order  to  form  the  Hessian,  the  second  partial  derivative  with  respect  to  p  of  Eq.  6.23 
must  be  calculated.  The  process  of  calculating  the  Hessian  begins  by  first  calculating 
the  first  partial  derivative  of  J  with  respect  to  p.  The  result  is  seen  in  Eq.  6.26. 


dJ 

dp 


EEE 

j= 1  k= 1  Z=1 


ei,k{^fl)Wj,k{uh) 


dp 


(6.26) 


Note  that  the  real  frequency  response,  Gj^{uo forest  is  not  affected  by  varying  model 
parameters.  Thus,  the  partial  derivative  of  the  frequency  response  error  £j,k  in  Eq.  6.26 
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can  be  represented  simply  as 


dejAluJfi)  _  /g  27) 

dp  dp 

It  is  important  to  define  that  dGl'^ is  the  jth  row  and  kth  column  of  the  frequency 
response  sensitivity  matrix  in  Eq.  6.29. 

d^rL  =  |  -  AylB  +  D)  (6-28) 

f)  A  F)¥t 

=  C(tu)fI  -A)-1  —  {iUfI  -  A)-lB  +  C{iujfI  -  H)”1  —  (6.29) 


Now,  in  light  of  Eq.  6.27,  the  first  partial  derivative  of  J  as  represented  in 
Eq.  6.26  can  be  represented  as 


dJ_ 

dp 


EEE 

j=  1  k= 1  1=1 


Ej,k{^fi )  (u Jfl ) 


dp 


(6.30) 


Taking  the  partial  of  Eq.  6.30  with  respect  to  p,  the  Hessian  Hj ^  is  formed,  as  seen 
in  Eq.  6.31. 


Hj:k  — 


i=i 


dp 


(6.31) 


Having  now  defined  Eq.  6.31,  the  definition  of  the  Cramer-Rao  bound  originally  stated 
in  Eq.  6.25  is  complete. 

Once  again,  the  preceding  discussion  is  strictly  applicable  only  to  an  LTI  system, 
and  is  not  applicable  to  LTP  models  such  as  a  helicopter  in  forward  flight.  Thus,  any 
parameter  identification  for  the  purpose  of  developing  a  controllable  linear  model  for 
such  an  LTP  system  will  not  be  verifiable.  Furthermore,  any  controller  based  on 
a  model  developed  using  those  parameters  may  therefore  be  inaccurate,  since  the 
parameters  may  not  accurately  represent  the  system  dynamics. 
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6. 6  Cramer-Rao  Methodology  for  LTP  Systems 

The  previous  Section  described  the  Cramer-Rao  methodology  for  an  LTI  system. 
However,  this  approach  is  not  applicable  to  LTP  systems,  such  as  a  helicopter  rotor 
in  forward  flight,  which  exhibits  periodic  behavior.  This  is  because  the  state  space 
representation  of  the  Cramer-Rao  bound  has  only  been  derived  for  an  LTI  system. 
For  this  approach  to  work  for  a  LTP  system  a  new  approach  must  be  created,  as  this 
has  not  been  done  to  date. 

To  derive  the  Cramer-Rao  bound  for  a  LTP  system  in  state  space  the  following 
procedure  will  be  taken.  First,  following  the  approach  in  Section  6.5,  the  linear  system 
must  be  described  in  state  space  form.  Wereley  [42]  has  developed  the  continuous  time 
Harmonic  Balance  State  Space  Model  (HBSSM),  which  provides  the  ability  to  analyze 
linear,  time  periodic  system  characteristics  using  techniques  developed  originally  for 
linear  time  periodic  systems.  Then,  following  a  brief  description  of  the  Harmonic 
Balance  State  Space  Model,  the  Cramer-Rao  bound  for  an  LTP  system  will  be  derived 
in  a  manner  similar  to  the  state  space  LTI  approach. 

6.6.1  LTP  State  Space  Operator.  A  linear,  time  periodic  system  is  said  to 
be  T-periodic  because  it  oscillates  with  period  T.  LTP  systems  are  T-periodic  due 
to  either  system  dynamics  or  an  input  signal  modulated  at  up,  which  is  the  system 
pumping  frequency ,  or  fundamental  frequency.  The  pumping  frequency  defines  the 
amplitude  and  period  at  which  the  system  modulation  takes  place.  This  concept  is 
illustrated  in  Figure  6.2,  where  a  simple  LTP  system  is  represented  as  a  time  varying 
gain,  1  —  2 (3  cosujpt  that  modulates  the  output  at  the  pumping  frequency  up. 

The  plant  dynamics  matrix,  Aft),  being  T-periodic  has  the  property 

Aft  +  T)  =  Aft)  Vi  G  (— oo,  oo)  .  (6.32) 

If  the  portion  of  the  matrix  defined  over  the  fundamental  interval  is  only  considered, 
then  Aft)  is  considered  bounded  such  that  A  ft)  €  L2[0,T\.  Considering  this  property, 
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1— 2cos  fiwpi 


u(t) 


y(t) 


Figure  6.2:  Example  of  a  simple  LTP  system  represented  by  a  modulating  gain. 

[42] 

the  state  space  model  of  a  T-periodic  system  is  described  as 


x(t )  =  A(t)x(t)  +  B(t)u(t ) 
y(t )  =  C(t)x(t )  +  D(t)u(t ) 


(6.33) 


noting  the  matrices  .£?(■),  C(-),  D(-)  are  also  T-periodic.  Now,  the  LTP  system 
represented  by  Eq.  6.33  can  be  excited  using  an  exponentially  modulated  periodic 
(EMP)  signal  of  the  form 


u(t)  =  umeSrnt  t  >  0 


(6.34) 


m=— oo 


noting  that  sm  —  s  +  jmup.  The  steady  state  response  to  an  EMP  signal  is 
represented  by  the  complex  Fourier  series 


m=— oo 


x(t)  = 

r 

oo 

Ht)  = 


Xm.& 


-Smt 


Sm.XmC 


(6.35) 

(6.36) 


m=— oo 


The  output  signal,  y(t),  is  represented  by  a  similar  expansion,  as  seen  in  Eq.  6.36. 
Likewise,  the  T-periodic  state  space  matrices  of  Eq.  6.33  can  be  represented  as  a 
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complex  Fourier  series,  as  seen  in  Eq.  6.37. 


OO 

A(t)  =  Y  AmejmuJpt 

m=— oo 


(6.37) 


where  u )p  is  the  LTP  system  pumping  frequency ,  or  fundamental  frequency.  The  matrix 
expansion  shown  in  Eq.  6.37  is  applied  identically  for  B(t),  C(t),  and  D(t )  of  Eq.  6.33. 
Now,  rewriting  the  state  space  representation  in  Eq.  6.33  with  Eqs.  6.36  and  6.37, 
results  in  the  Fourier  series  LTP  system  seen  in  Eq.  6.38. 

OO  OO  OO  OO  OO 

Y  snxneSnt  =  Y  A-ejnUpt  E  x™eSmt+  Y  B^™vt  E  u™eSmt 

n=— oo  m=—o o  m=— oo  m=— oo  m=— oo 

oo  oo 

E  AnXmeSn+mt  +  Y  BnUmeSn+rnt 

n,m=— oo  n,m=— oo 

oo  oo 

=  Y  An-mxmeSrit  +  Y  Bn_mumeSnt  (6.38) 

n,m=— oo  n,m=— oo 

Note  once  again  that  sn  =  s  +  jmop  and  sm  are  represented  likewise.  Using  the 
same  technique,  the  measurement  equation,  y(t )  of  Eq.  6.33  can  be  represented  in  a 
manner  similar  to  that  of  Eq.  6.38,  as  seen  in  Eq.  6.39. 

OO  OO  OO 

Y  yneSnt  =  Y  Cn-mxmeSnt+  Y  Dn-mumeSnt  (6.39) 

n=— oo  n,m=—o o  n,m=— oo 


The  principle  of  harmonic  balance  [12,42]  can  be  applied  to  Eqs.  6.38  and  6.39 
by  moving  all  terms  on  the  right-hand  side  of  the  equations,  and  then  multiplying 
through  both  equations  by  e~st .  The  resulting  equations  generate  the  modified  state 
space  representation  of  the  LTP  system  as  seen  in  Eq.  6.40. 


OO 


OO 


^  ^  ^4 n—m^m  “1“  ^  ^  ^n—m 

m=— oo  m=— oo 

oo  oo 

Vn  ^  ^  Cn—m%m  4“  ^  ^  D 


n—m^m  \  /  J 

m=— oo  m=— oo 


(6.40) 
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The  harmonic  state  space  system  in  Eq.  6.40  can  be  simplified  into  the  Toeplitz  form, 
shown  in  Eq.  6.41.  This  result  is  defined  as  the  Harmonic  Balance  State  Space  Model 
(HBSSM). 


sx  =  {A  —  M)x  +  Bu 

y  =  Cx  +  Vu  (6.41) 


where  x,  u ,  and  y  represent  the  doubly  infinite  state,  control,  and  output  vectors, 
respectively,  defined  in  terms  of  modulated  complex  Fourier  series  coefficients,  as 
shown  in  Eq.  6.42. 


X_2 

'U-2 

y —2 

X-l 

u- 1 

y  i 

x0 

,  u  = 

u0 

,  y  = 

2/o 

Xi 

Ul 

2/1 

x2 

u2 

2/2 

(6.42) 


The  Toeplitz  form  of  each  of  the  state  space  matrices  A,  b,  c,  and  v  represented  in 
Eq.  6.41  is  expressed  as  a  doubly  infinite  block  of  the  complex  Fourier  coefficients, 
An\n  G  Z.  Eq.  6.43  represents  the  T-periodic  dynamics  matrix,  A(t),  of  Eq.  6.40  in 
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Toeplitz  form. 


Ao 

1 

to 

1 

CO 

A— 4 

A\ 

Ao 

1 

to 

1 

CO 

a2 

Ai 

Ao 

1 

to 

A3 

a2 

A ! 

Ao 

A_! 

A4 

A3 

A2 

A, 

Aq 

(6.43) 


The  matrix  M  in  Eq.  6.41  is  a  doubly  infinite  block  diagonal  containing  multiples  of 
the  pumping  frequency,  which  is  defined  as 

J\f  =  blkdiag{jnujpI}  Vn  <E  Z  .  (6.44) 

Using  the  definition  of  the  Toeplitz  form  of  the  harmonic  state  space  model  in  Eq.  6.41, 
the  linear  operator  ,g,  which  relates  input  to  output,  can  be  defined  as 


y  =  G(s)u 


(6.45) 


where 


G{s)  =  C[sZ-(A-N)]-lB  +  V  .  (6.46) 

6. 6.2  Cramer-Rao  Definition  for  LTP  Systems.  In  this  Section,  the  Cramer- 
Rao  methodology  analogous  to  that  described  in  Section  6.5  will  be  developed  for  a 
LTP  system  using  the  Harmonic  Balance  State  Space  Model,  Eq.  6.41.  This  method¬ 
ology  provides  a  convenient  method  to  quickly  construct  the  bound  for  the  individual 
model  parameters,  but  has  been  restricted  to  LTI  systems.  This  restriction  is  primar¬ 
ily  due  to  the  inability  to  accurately  describe  a  LTP  system  as  a  true  linear  operator 
having  input  and  output  signal  spaces  of  the  same  dimension.  To  accurately  under- 
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stand  how  the  linear  operator  theoretic  description  in  Section  6.5  is  applied  to  a  LTP 
system,  it  is  important  to  first  describe  the  differences  in  the  signal  spaces  of  LTI  vs. 
LTP  systems. 

A  LTI  system  having  a  sinusoidal  input  signal  of  frequency  Uf  will  map  an 
output  signal  of  un  but  with  different  magnitude  and  phase  than  the  input  signal. 
Thus,  the  signal  spaces  are  identical,  as  the  input  and  output  are  both  complex 
exponential  sinusoids  of  frequency  a if.  This  response  for  a  LTI  system  is  clearly  seen 
in  a  Bode  plot,  as  phase  and  magnitude  of  the  system  are  plotted  over  the  entire 
range  of  input  frequencies  ,  u)f  —  0  . . .  oo. 

A  LTP  system  differs  from  a  LTI  system  in  that  the  system  parameters  vary 
periodically  rather  than  remain  constant,  or  time  invariant.  Consider  a  LTP  system 
defined  as  such  by  having  a  time  periodic  input  matrix,  B.  In  this  case  the  steady  state 
forced  response  signal  of  a  LTP  system  will  be  composed  of  the  sinusoidal  input  signal 
u>f  and  the  pumping  frequency  input  signal  up,  as  seen  in  Figure  6.2.  In  this  Figure 
an  input  signal,  u(t),  is  convoluted  with  a  time  periodic  gain  1  —  2j3cosu>pt.  Thus, 
the  output  signal,  y(t),  is  a  convolution  of  signals  best  described  by  the  Fourier  series 
containing  the  frequencies  Uf  ±  nu>p\ujf  G  C ,  n  G  Z,  which  is  an  infinite  dimension 
signal  space.  This  is  depicted  in  Figure  6.3.  Clearly,  the  input  and  output  spaces 
are  not  the  same  dimension  as  in  the  case  of  the  LTI  system  spaces.  This  difference 
in  signal  spaces  differentiates  the  LTP  system  from  LTI  systems,  as  the  mapping  of 
signal  spaces  from  input  to  output  is  not  an  isomorphism ,  being  one-to-one  and  onto, 
but  rather  a  one-to-many  mapping  .  It  is  for  this  reason  that  in  this  form,  the  LTP 
system  is  not  considered  a  true  linear  operator. 

Wereley  [42]  developed  the  concepts  of  Exponentially  modulated  periodic  signal, 
as  seen  in  Equation  6.34,  to  rectify  the  signal  space  imbalance  in  LTP  systems.  The 
EMP  signal  created  an  input  space  of  equal  dimension  to  the  output  signal  space, 
which  formed  the  foundation  for  the  LTP  operator  described  in  section  6.6.1.  By  using 
the  LTP  harmonic  state  space  model  from  definition  6.45  the  input-output  relationship 
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Linear  Time-Constant  System: 

Sinusoidal  Input  Sinusoidal  Response 


Figure  6.3:  Multiharmonic  Response  of  an  LTP  System  [20]. 
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between  the  input  harmonics,  un\n  E  Z,  and  the  output  harmonics,  yn\n  e  Z.  Thus, 
a  true  linear  operator  for  LTP  systems  is  defined. 

As  the  HBSSM  is  a  true  linear  operator,  the  Cramer-Rao  lower  bound  can  be 
described  using  the  state  space  matrices  A ,  B ,  C,  and  V  from  Eq.  6.41  in  a  modified 
form  of  Equation  6.29,  as  seen  below. 


dQ(tujf) 

dp 


^  [C(tufI  -  AUU)-XB  +V] 

BA 

C(lu}I  -AHATy'—iLujfl  -  AHN)~lB 

BB 

+  C{uujfI  -  AHN)-1  — 


where  the  matrix  AHAf  is  defined  as 


(6.47) 


(6.48) 


AHN  =  A  -  N 


(6.49) 


The  above  Equations  6.48  and  6.49  redefine  the  partial  derivative  of  the  response  error 
from  Equation  6.27  to  an  expression  containing  a  doubly  infinite  block  of  complex 
Fourier  coefficients  in  Toeplitz  form.  To  complete  the  computation  of  the  Cramer- 
Rao  lower  bound  for  LTP  systems  defined  in  block  Toeplitz  form,  Equation  6.48  is 
used  to  in  Equation  6.50  to  define  the  Hessian,  H,  in  a  manner  analogous  to  Equation 
6.31. 


H-j,k  —  kEj,fc(a>;) 


i=i 


BQ j,k{LU  ft) 
Bp 


(6.50) 


Having  defined  the  Hessian  for  the  LTP  system  in  block  Toeplitz  form  in  Equation 
6.50,  the  Cramer-Rao  lower  bound, cn,  is  now  defined  for  LTP  systems  as 


ClZj  ~ 


(6.51) 
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u(t) 


V 


Figure  6.4:  Example  LTP  system  with  time  invariant  A,C  matrices  and  time  peri¬ 
odic  B(t)  matrix. 

6.7  Illustrative  Example 

Having  defined  the  Cramer-Rao  lower  bound  for  LTP  systems  in  Equation  6.51, 
an  example  will  now  show  how  the  method  can  be  applied  to  LTP  systems  to  validate 
identified  system  parameters.  For  this  example,  the  simple  LTP  system  described  in 
Figure  6.4  will  be  used  to  evaluate  the  validity  of  the  system  parameter  (. 

For  the  example,  the  system  depicted  in  Figure  6.4  is  considered  with  time 
invariant  plant  and  output  matrices  Aa  and  CQ  and  a  T- periodic  control  matrix  B(t). 
The  system  matrices  are  defined  as  seen  in  Equation  in  6.52 


0  1 

0 

B{t)  = 

a  =  (h) 

1  —a 

^ n  — 

1  —  2  /3cosuipt 

\  11  / 

(6.52) 


where  the  values  of  the  system  are  defined  as  ujv  =  1 ,  un  =  .5,  (  =  .3,  a  =  — 1, 
and  /3  —  1.  As  the  B(t)  matrix  is  to  be  described  by  way  of  the  Fourier  series  for 
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use  in  the  Toeplitz  transform,  it  is  therefore  expressed  as  in  Equation  6.53. 


B(t)  = 


,  0  ,  B-i  ,  Bq  ,  B\  ,  0  , 


,  0  , 


0 

-p 


j  Bq  , 


0 

-p 


,  0  , 


(6.53) 


These  matrices  are  used  to  create  the  ffBSSM  of  this  model  in  the  manner  described 
in  Section  6.6.1.  It  is  important  to  note  that  while  the  number  of  harmonics  of  the 
HBSSM  is  doubly  infinite,  for  use  in  a  digital  computer,  the  system  is  truncated  to 
include  only  n  positive  harmonics  with  an  equal  number  of  negative  harmonics  along 
with  the  the  zero  harmonic.  To  ensure  accurate  results  in  this  example,  n=10  positive 
harmonics  are  used.  The  rationale  for  this  selection  is  discussed  below. 

The  Toeplitz  block  form  of  Fourier  coefficients  is  key  to  developing  a  state  space 
form  of  the  system  matrices  that  define  a  model.  These  matrices,  as  described  previ¬ 
ously,  are  considered  doubly  infinite  as  they  contain  a  doubly  infinite  representation 
of  the  Fourier  series.  This,  however  is  untenable  for  application  in  code  development 
on  a  digital  computer.  For  code  development  to  take  place,  the  representative  Fourier 
series  of  the  Toeplitz  transform  must  be  truncated  to  contain  N  harmonics.  To  de¬ 
termine  the  value  of  N,  the  number  of  harmonics  was  increased  in  the  system  as  to 
identify  convergence  in  the  compared  system  responses.  Based  on  the  experience  of 
Hwang  [12],  n=10  was  determined  to  produce  adequate  convergence  in  system  output 
response. 

By  creating  the  HBSSM  for  this  example,  the  block  Toeplitz  matrices  a,  b,  c,  and  v 
along  with  the  partial  derivatives  and  are  then  used  in  Equation  6.48  to  generate 
the  Hessian  of  Equation  6.50.  The  weighting  matrix,  W  in  Equation  6.50  is  formed 
by  using  an  input  noise  spectral  density,  Sv  —  1,  for  the  input  noise  covariance  ma¬ 
trix,  R.  Thus,  for  this  example,  the  weighting  matrix  is  unitary  as  W  =  R =  1  . 

As  this  example  is  determining  the  Cramer-Rao  bound  for  the  system  parameter  £, 
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CRAMER-RAO  PLOT 


Figure  6.5:  Cramer-Rao  bound  of  the  system  in  Equation  6.52  with  respect  to  the 
parameter  (. 

the  partial  derivatives  are  taken  with  respect  to  (.  After  determining  the  Hessian, 
the  Cramer-Rao  bound  for  (  is  determined  by  use  of  Equation  6.51.  A  plot  of  the 
Cramer-Rao  bound  for  (  is  presented  over  a  frequency  range  toj  G  (— qf  >  qf  ]  ■  The 
Cramer-Rao  bound  presented  in  Figure  6.5  is  based  solely  on  the  first  system  har¬ 
monic,  which  is  the  dominant  mode  of  the  system.  The  results  of  this  plot  indicate 
that  the  lowest  Cramer-Rao  bound  about  the  parameter,  (  is  a  standard  deviation 
of  ±  0.004,  representing  a  total  bound  of  0.008.  Additionally,  this  value  occurs  at  an 
input  frequency  of  Uf  =  0.5.  In  terms  of  validating  the  test  parameter  by  way  of  the 
Cramer-Rao  bound,  any  value  of  (  derived  empirically  that  has  a  Cramer-Rao  bound 
significantly  larger  than  the  smallest  Cramer-Rao  bound  may  be  invalid  and  warrants 
further  investigation. 


85 


Identified  Values  of  Parameter Zefa  with  Cramer-Rao  Bounds 
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Figure  6.6:  Values  of  (  derived  by  system  identification  at  input  frequency  ujf  with 
superimposed  Cramer-Rao  bounds. 

An  example  of  using  the  Cramer-Rao  bound  to  validate  parameter  estimates 
is  presented  in  Figure  6.6.  Here  parameter  estimates  of  (  are  presented  at  select 
input  frequencies,  Uf  =  0,0.05,0.1,0.2,0.3,0.4,0.5  rad/s.  As  a  note,  the  parameter 
estimates  of  (  were  performed  by  using  a  frequency  domain  parameter  identification 
method  for  LTP  systems  developed  by  Hwang  [12],  The  output  data  used  for  these 
parameter  estimates  was  corrupted  by  the  zero  mean  Gaussian  white  measurement 
noise,  u,  of  Figure  6.4  having  spectral  density  Sv  =  1.  This  spectral  density  for  the 
input  noise  was  chosen  to  match  the  value  used  to  define  the  weighting  matrix  W 
of  the  Cramer-Rao  bound.  The  parameter  estimates  of  (  are  presented  with  their 
corresponding  Cramer-Rao  bound  from  the  data  presented  in  Figure  6.5.  The  values 
of  (  are  presented  at  each  frequency  used  for  system  identification  to  compare  their 
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accuracy  to  the  true  value  of  (  =  0.3  in  addition  to  showing  the  size  of  the  corre¬ 
sponding  Cramer-Rao  bound.  As  stated  in  Section  6.5,  a  large  Cramer-Rao  bound 
indicates  poor  parameter  estimation  performance,  and  thus  by  reviewing  Figure  6.6 
one  can  identify  accurately  and  inaccurately  identified  values  of  (.  Thus,  the  values 
of  (  in  Figure  6.6  that  fall  in  the  input  frequency  range  of  0.3  <  ujf  <  0.5  have  both 
identified  values  of  (  near  to  the  true  value  of  0.3,  and  demonstrate  small  Cramer-Rao 
bounds,  indicating  a  small  standard  deviation  of  the  identified  value.  These  values 
are  therefore  accurate  and  should  be  considered  for  use  in  the  parametric  system 
model.  Likewise,  the  identified  values  of  (  that  fall  in  the  input  range  of  0  <  Uf  <  0.3 
have  excessively  large  bounds,  indicating  poor  identification  of  which  is  evident  by 
their  values  straying  from  the  true  value  of  (  =  0.3.  The  data  in  this  range  should 
therefore  not  be  used.  If  the  bound  of  this  parameter  grows,  it  can  indicate  poor 
data  correlation  and  thus  warrants  investigation  of  the  test  data  used  for  parameter 
identification.  Finally,  it  is  important  to  note  that  this  plot  presents  not  only  the  best 
identified  value  of  (  based  on  the  Cramer-Rao  bound,  but  also  the  value  of  the  input 
frequency,  Uf  at  which  this  occurs.  This  can  be  used  to  specify  the  ideal  input  fre¬ 
quency  range  to  use  when  performing  system  identification  that  will  ensure  accurate 
parameter  results. 

6.8  Concluding  Remarks 

A  derivation  of  the  theory  and  methodology  required  to  generate  the  Cramer- 
Rao  lower  bound  for  a  specified  parameter  in  a  linear,  time  periodic  (LTP)  system 
in  state  space  form  has  been  presented.  This  development  now  makes  it  possible  to 
determine  the  bounded  standard  deviation  of  a  system  parameter  which  has  been 
estimated  using  any  system  identification  technique.  The  Cramer-Rao  lower  bound 
represents  the  standard  deviation  based  on  using  an  optimal  estimator,  thus  providing 
a  true  measure  of  the  accuracy  of  the  estimate. 

Through  an  illustrative  example,  it  has  been  shown  that  in  a  LTP  system,  the 
frequency  at  which  the  parameter  estimation  is  performed  is  critical  to  the  confidence 
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that  one  can  attribute  to  the  estimate.  Therefore,  knowledge  of  the  Cramer-Rao  lower 
bound  provides  those  performing  the  system  identification  with  valuable  insights  on 
how  best  to  obtain  estimates  of  system  parameters.  The  importance  of  this  knowledge 
has  also  been  shown  to  be  particularly  important  when  there  is  noise  in  the  data  used 
to  perform  the  parameter  estimates. 

As  noted  in  the  Introduction,  minimizing  the  one-per-rev  vibrations  of  helicopter 
rotors  is  important  to  the  safe,  efficient  operation  of  helicopters.  The  development  of  a 
methodology  for  determining  the  Cramer-Rao  lower  bound  for  LTP  systems  provides 
a  means  for  assessing  the  quality  of  parameterized  models  developed  for  rotor  smooth¬ 
ing.  Further  research  along  this  line  will  focus  on  using  the  Cramer-Rao  methodology 
to  improve  the  efficiency  of  rotor  smoothing  methods.  One  possibility  for  reducing 
the  number  of  flights  and  manual  adjustments  is  to  introduce  actuators  in  the  rotor 
system.  These  actuators  would  be  controlled  by  an  on-line  system  which  performs 
continuous  system  identification,  thereby  providing  verifiably  accurate  control.  Such 
a  system  would  alleviate  the  repetitiveness  of  the  current  rotor  vibration  adjustment 
process  by  providing  continuous  adjustments  for  vibration  reduction  during  flight. 


VII.  Rotor  Vibration  Smoothing  Using  Cramer-Rao 

Parameter  Validation 


The  overall  intent  of  this  research,  as  stated  in  chapter  II,  is  to  develop  an  op¬ 
timal  rotor  smoothing  approach  to  reduce  out  of  plane  vibrations  generated  by 
asymmetrical  aerodynamic  lift  of  individual  blades.  The  proposed  rotor  smoothing 
method  is  based  upon  an  framework  containing  the  following  three  tenets: 

1.  Perform  system  identification  to  populate  the  rotor  system  parametric  model. 

2.  Validate  the  accuracy  of  the  identified  parameters  of  the  rotor  system  model. 

3.  Produce  a  vibration  control  solution  using  linear  optimal  methods  based  upon 
the  validated  rotor  parameters  to  reduce  out  of  plane  vibrations  to  an  acceptable 

level. 

The  previous  chapters  have  provided  the  necessary  tools  to  accomplish  this 
approach,  with  an  emphasis  on  the  Cramer-Rao  bound  LTP  systems  to  provide  system 
parameter  validation.  This  chapter  will  culminate  the  works  of  the  previous  chapter  to 
produce  a  LTP  Linear  Quadratic  Regulator  (LQR)  based  upon  Cramer-Rao  validated 
system  parameters.  Through  this  approach,  the  utility  of  the  Cramer-Rao  bound  will 
be  evident,  as  it  can  detect  poorly  identified  system  parameters,  which  ultimately 
lead  to  a  poorly  performing  LQR  vibration  controller. 

7.1  Needed  Improvements  in  Rotor  Smoothing  Algorithms 

Rotor  smoothing  algorithms  for  rotors  are  not  in  short  supply,  as  outlined  in 
chapter  II.  However  an  accurate  approach  is  still  elusive.  This  is  due  in  large  part  to 
a  smoothing  approach  based  upon  a  flawed  internal  model  of  the  rotor  system.  The 
primary  flaw  in  these  rotor  smoothing  approaches  is  the  use  of  a  linearized  model 
that  does  not  capture  the  periodic  nature  of  a  rotor  system  in  forward  flight.  While 
it  is  evident  that  a  linearized  flapping  rotor  blade  model  contains  periodic  terms  in 
the  plant  and  input  matrices,  A(t)  and  B(t)  respectively,  the  current  rotor  smoothing 
approaches  using  linear  models  average  the  periodic  terms  in  order  to  produce  a  LTI 


system.  This  approach  has  been  done  in  order  to  use  existing  control,  system  identi¬ 
fication,  and  parameter  validation  techniques  that  are  valid  for  LTI  systems.  While 
this  simplification  has  allowed  for  the  development  of  rotor  smoothing  algorithms,  it 
has  incurred  the  inaccuracies  associated  with  using  a  system  that  does  not  capture 
the  time  periodic  nature  of  the  flapping  blade  in  forward  flight. 

An  additional  flaw  in  existing  rotor  smoothing  approaches  is  the  lack  of  vali¬ 
dated  parameters  in  the  internal  rotor  system  model  used  to  develop  the  smoothing 
solution.  Rotor  smoothing  algorithms  that  use  an  internal  system  model  to  develop 
the  vibrational  smoothing  solution  rely  on  the  measured  response  of  the  actual  system. 
This  data  is  in  turn  used  to  construct  a  representative  internal  model  used  to  derive 
the  rotor  smoothing  solution,  as  seen  in  chapter  II.  The  representative  model  can  be 
either  a  parametric  model,  usually  in  State  Space  representation,  or  a  non-parametric 
model  such  as  a  frequency  response  function.  For  a  parametric  model,  the  validity 
of  the  individual  parameters  can  be  evaluated  for  accuracy  before  use  using  existing 
approaches,  such  as  the  Cramer-Rao  approach  covered  in  detail  in  chapter  VI.  For 
non-parametric  rotor  smoothing  algorithms,  such  as  the  U.S.  Army’s  AVA,  a  para¬ 
metric  validation  approach  is  reduced  to  validating  the  single  parameter  representing 
the  frequency  response,  which  further  reduces  the  ability  to  detect  modeling  errors 
as  there  are  no  physical  parameters  that  can  be  accurately  validated.  Regardless  of 
modeling  approach,  parameter  validation  is  rarely  ever  used  in  rotor  smoothing  and 
thus  modeling  inaccuracies  that  could  otherwise  be  detected  by  a  parameter  valida¬ 
tion  approach  are  introduced  into  the  rotor  vibration  reduction  solution.  Maine  and 
Iliff  [34]  note  that  if  highly  accurate  estimates  cannot  or  are  not  distinguished  from 
worthless  estimates,  to  be  safe  all  estimates  must  be  considered  suspect  or  moreover, 
worthless.  On  that  note,  if  the  smoothing  algorithm  is  to  be  accurate,  a  model  based 
upon  validated  parameters  is  essential  to  a  well  performing  smoothing  algorithm. 

The  following  Section  will  describe  the  methodology  used  to  overcome  these 
deficiencies. 
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1.2  Outline  of  Cramer-Rao  Validated  Controller  Development 

The  work  of  this  chapter  is  to  describe  how  to  smooth  main  rotor  system  out 
of  plane  vibrations  through  a  LTP  LQR  controller.  In  addition  to  the  controller 
development,  this  chapter  will  demonstrate  the  importance  of  the  Cramer-Rao  lower 
bound  for  LTP  system  parameters  by  comparing  the  performance  of  several  LQR 
controllers.  Here,  each  controller  is  developed  based  on  upon  validated  parameters 
having  a  Cramer-Rao  bound  different  from  another.  For  the  controller  comparison, 
the  Cramer-Rao  bounds  will  be  increased  in  size  for  each  successive  controller.  As 
shown  in  chapter  VI,  the  validity  of  system  parameters  derived  from  a  parameter 
identification  scheme  is  garnered  from  the  magnitude  of  the  Cramer-Rao  bound.  This 
is  due  to  the  fact  that  the  Cramer-Rao  bound  depicts  the  standard  deviation  of  the 
identified  parameter.  As  the  performance  of  the  LQR  is  directly  related  to  the  model 
perturbations,  which  arise  from  identified  system  parameters  having  varying  levels  of 
accuracy,  a  method  to  quantify  the  level  of  parameter  perturbations  would  allow  one 
to  determine  the  controller  effectiveness  simply  by  reviewing  the  magnitude  of  the 
bound  itself.  Thus  ,  the  Cramer-Rao  bound  will  demonstrate  that  poorly  identified 
system  parameters  will  lead  to  unsuitable  controller  demands  on  control  inputs. 

In  order  to  achieve  the  above  stated  goals  the  following  steps  will  be  accom¬ 
plished  in  successive  order  in  this  chapter: 

1.  Describe  the  rotor  system  LTP  equations  of  motion  based  upon  a  rigid  blade 
model,  with  4  blades  in  total. 

2.  Develop  the  LTP  LQR  controller  for  out  of  plane  vibration  reduction 

(a)  Describe  the  properties  of  a  LTP  LQR  controller 

(b)  Outline  the  tracking  regulator  design  used  to  reduce  out  of  plane  vibrations 

(c)  Explain  selection  of  time  periodic  gain  harmonics 

3.  Outline  the  computation  of  each  controller  based  upon  increasing  Cramer-Rao 
bounds. 
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4.  Compare  the  performance  of  each  controller  developed  in  the  above  step  to 

demonstrate  the  effects  of  large  Cramer- Rao  bounds  on  controller  performance 

The  following  Sections  will  address  the  individual  goals  in  succession,  beginning  with 
the  derivation  of  the  rotor  system  equations  of  motion. 

7.3  Rotor  System  LTP  Equations  of  Motion 

As  stated  in  the  introduction  to  this  chapter,  each  rotor  smoothing  algorithm 
relies  upon  an  internal  rotor  model  in  order  to  compute  the  rotor  smoothing  solution. 
For  this  work  a  simplified  four  bladed  rotor  system  is  used  to  represent  the  dynamics 
of  the  flapping  rotor  blades  in  forward  flight.  As  this  work  is  only  considering  rotor 
smoothing  and  not  higher  harmonic  control,  rigid  blades  are  adequate  and  therefore 
will  be  used  to  represent  the  individual  blades  in  the  rotor  system  model.  Furthermore, 
for  the  sake  of  simplification,  each  blade  will  be  considered  without  a  flap  hinge  offset 
or  spring  restraint. 

The  rotor  system  will  be  derived  by  first  developing  the  equations  of  motion  of 
an  individual  blade  having  the  characteristics  described  in  the  last  paragraph.  Once 
the  individual  blade  model  is  derived  a  system  comprised  of  four  identical  blades  will 
be  assembled,  with  each  blade  phased  90  degrees  apart,  as  in  the  configuration  of  the 
AH-64  as  seen  in  Figure  7.1  or  similar  4  bladed  rotor.  The  derivation  for  a  single  rigid 
blade  that  follows  is  derived  from  the  work  of  Johnson  [16],  and  is  presented  below. 

The  rigid  blade  flapping  equations  of  motion  are  now  considered  by  first  taking 
the  equilibrium  of  the  inertial  and  aerodynamic  moments  about  the  flapping  hinge,  as 
seen  in  Figure  7.2.  Considering  the  blade  mass  element,  m  dr ,  at  the  radial  position 
r  the  the  three  sectional  forces  acting  on  the  mass  are  described  as  follows: 

1.  The  Inertial  force,  described  as  mz  =  mr/3 

2.  The  Centrifugal  force,  described  as  mfl2r 

3.  The  Aerodynamic  force,  Fz 
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Figure  7.1:  Diagram  of  AH-64  rotor  system  [14]. 
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where  z  =  (3r  describes  the  out  of  plane  deflection  of  the  flapping  rotor  blade.  As 
each  of  these  forces  are  considered  in  terms  of  moments  about  the  flapping  hinge,  the 
radial  distances  for  each  force  are  considered  as 

1.  The  inertial  force  has  a  moment  arm  r  about  the  hinge 

2.  The  centrifugal  force  acts  radially  outward  form  the  blade  with  moment  arm 

z  =  r/3 

3.  The  Aerodynamic  force  Fz  has  a  moment  arm  r  about  the  hinge 

The  three  forces  act  to  balance  the  blade  in  equilibrium  when  the  system  rotates  with 
speed  fi.  Here,  the  aerodynamic  force  Fz  is  considered  as  the  lift  of  the  individual 
section  of  the  blade,  initiating  the  upward  flap  motion  registered  in  terms  of  the  flap 
angle  (  /3  ).  The  centrifugal  force  mfl2r  and  inertial  force  mrf3  act  to  oppose  the 
flapping  motion  caused  by  the  aerodynamic  force  at  each  blade  section. 

The  equilibrium  condition  is  generated  taking  the  sum  of  the  moments.  The 
moments  in  this  case  will  be  equal  to  zero  as  there  is  no  blade  hinge  spring  considered 
in  this  case.  The  moments  are  generating  by  integrating  the  sectional  forces  over 
the  entire  blade  span  from  root  to  tip  while  taking  the  product  with  respect  to  the 
corresponding  moment  arm  at  the  location  of  the  sectional  force.  This  operation  is 
represented  in  Equation  7.1. 

pR  pR  pR 

/  mrf3r  dr  +  mfl2r(r/3)  dr  —  /  Fzr  dr  =  0  (7-1) 

Jo  Jo  Jo 

The  representation  of  the  moment  equilibrium  in  Equation  7.1  can  be  adjusted  to  put 
the  aerodynamic  moments  on  the  right  hand  side  of  the  equation,  as  seen  in  Equation 

pR  pR 

/  mr/3r  dr  +  mtt2r(r(3 )  dr  =  fQ  Fzr  dr  (7-2) 

Jo  Jo 
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Figure  7.2:  Rigid  Rotor  Blade  Flapping  Moments  [16]. 
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7.2.  The  equilibrium  equations  can  be  further  simplified  by  defining  the  blade  moment 
of  inertia  about  the  flapping  hinge  (Ifl)  as  represented  in  Equation  7.3. 


rR 

h  —  I  r2m  dr  (7.3) 

Jo 


Thus,  Equation  7.2  reduces  to  Equation  7.4  by  way  of  Equation  7.3,  as  seen  below. 


fj  +  P 


(7.4) 


As  a  note,  the  above  equation  is  transformed  to  use  dimensionless  quantities  for  12 
and  R.  Additionally,  the  air  density  coefficient  (  p  )is  normalised  for  use  in  the  Lock 
number (7  derivation.  The  Lock  number  represents  the  ratio  of  aerodynamic  forces  to 
inertial  forces  in  an  dimensionless  parameter.  It  is  important  to  note  that  the  Lock 
number  contains  the  sole  influence  on  flap  motion  by  way  of  the  air  density,  as  seen 
in  the  Lock  number  definition  below.  This  will  play  a  major  part  in  defining  the  LTP 
equations  for  flapping  blade. 

pacR 4 


As  a  note,  the  parameters  a  and  c  in  the  Lock  number  equation  above  represent 
the  blade  section  two  dimensional  lift  curve  slope(a)  and  the  blade  chord  widthic) 
respectively.  The  Lock  number  equation  can  be  combined  with  Equation  7.4  to  derive 
the  second  order  flapping  blade  equation,  as  seen  in  7.6. 

f1  F 

(3  +  (3  =  7  /  r—dr  =  ^MF  (7.6) 

Jo  ac 


In  order  to  define  the  total  blade  flapping  equation  the  right  hand  side  of  Equa¬ 
tion  7.6  must  be  defined.  The  right  hand  side  of  this  equation  is  termed  the  aero¬ 
dynamic  flap  moment (  MF  ),  which  is  simply  the  normalised  aerodynamic  force,  ^ 
acting  normal  to  the  blade  at  the  radial  position,  r.  The  normalised  aerodynamic 
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force  can  be  more  accurately  derived  in  terms  of  the  tangential  and  perpendicular  air 
velocities  at  the  blade  segment,  (  Up  )  and  (  Up  )  respectively,  as  seen  in  Equation 
7.7. 


—  =  -  =  Uufo  -  uPuT ) 

ac  ac  2 


(7.7) 


where  the  blade  pitch  angle  is  denoted  by  (  0  ).  By  using  the  normalised  aerodynamic 
force  of  Equation  7.7  the  aerodynamic  flap  moment  can  be  described  as 


MP  = 


f1  Fz 

/  r  —  dr 

I  o 


-I 


JO 

1 

=  r  — 
2 


ri(U}B  -  UPUT) 


(7.8) 

(7.9) 


(r  +  psinip)2Q  —  (A  +  r/3  +  n  (3  cosi/j)(r  +  psimp)  dr{ 7.10) 


The  term  ip  in  the  above  equation  refers  to  the  dimensionless  time  variable  for  the 
rotor  azimuth(  ip  ),  which  is  related  to  the  rotational  velocity  by  ip  =  Qt  .  Also, 
the  terms  p  and  A  refer  to  the  rotor  advance  ratio  (  p  )  and  rotor  inflow  ratio  (  A  ) 
respectively.  These  terms  are  used  to  describe  the  forward  speed  of  the  rotor  system 
and  the  ratio  of  the  tangential  air  inflow  due  to  forward  velocity,  V,  versus  the  inflow 
due  to  the  rotating  blades  of  the  rotor  system,  respectively.  As  a  note,  this  model 
assumes  linear  blade  twist  and  uniform  inflow.  Performing  the  integration  of  Equation 
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7.10  results  in  the  following  expression  for  the  flapping  moment 
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(7.11) 

(7.12) 


where  0con  is  the  control  input  for  blade  pitch (  Ocon  )  and  Qtw  is  the  blade  twist  ( 
(~)tw  ).  Therefore  ,  by  using  the  above  dehned  representation  of  the  flapping  moment, 
Mp,  Equation  7.6  can  be  revised  as 


P  +  P 
P  +  P 


'fMp 
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+  M\\  +  M +  Mpf3 


(7.13) 


By  reviewing  Equation  7.14  with  consideration  of  the  coefficients  making  up  the  flap¬ 
ping  moment  Mp  in  Equation  7.13  it  is  clear  that  the  second  order  flapping  equations 
of  motion  are  time  periodic.  This  is  due  to  the  sin  and  cosine  terms  in  the  individual 
coefficient  that  make  up  Mp.  Equation  7.14  can  be  further  simplified  by  linearising 
the  terms  of  Mp  and  combining  like  terms  to  result  in 


P+  |  f  1  +  ^hsin^  )  @  + 


rf  (  Zj. 

1  +  —  (  -  jJL  cos  xp  +  [i2  sin  2ip 


P  =  fW)  .  (7.14) 


It  is  noted  with  balancing  of  terms  that  the  terms  of  MF  that  are  left  on  the  right 
hand  side  define  a  periodic  forcing  function,  f{ip)  as  seen  in  Equation  7.15. 
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Further  simplifications  can  be  made  by  assuming  zero  blade  twist,  0^,  and  neglecting 
A  as  it  can  be  assumed  as  a  plant  disturbance  in  this  case. 

The  linearized  equations  of  motion  represented  in  Equation  7.14  can  be  trans¬ 
formed  into  state  space  form  represented  in  7.16  having  states  f3,  (3  and  having  control 
input  0con. 


y{i 0  =  ctV’MVO  (7.16) 


where  the  feedforward  matrix  D( ip)  is  omitted  as  there  is  110  feedforward  in  this 
system.  As  noted  above,  the  system  is  linear  time  periodic,  and  is  described  in  state 
space  form  in  term  of  the  time  periodic  matrices  A(i(j)  and  B(iJ>)  as 


A(jp)  = 

B(i 0  = 


0  1 
-  {1  +  Kl/xcos-^  +  p2  sin  2-0)}  — |  (l  +  |//sin -0) 

0 

|  (l  +  /j2  +  |p  sin  ^  —  /i2  cos  2^) 


(7.17) 

(7.18) 


The  output  of  this  system  is  only  the  state  /3,  which  is  represented  by  the  output 
matrix,  C( ij}) 


cw 


1  0 


(7.19) 
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In  order  to  transform  the  linear  time  periodic  system  in  to  harmonic  balance  state 
space  form,  the  Fourier  series  coefficients  of  the  time  periodic  matrices  in  Equations 
7.18  and  7.18  must  first  be  identified.  This  is  seen  in  Equations  7.20  and  7.21  ,  for  A 
and  B  respectively.  As  the  output  matrix  C  has  no  periodic  terms,  it  will  have  only 
the  zero  harmonic  term,  Co- 
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It  is  important  to  note  that  the  terms  denoted  by  an  asterix  in  Equations  7.20  and  7.21 
are  the  complex  conjugates  of  the  respective  matrix.  In  light  of  the  above  comments, 
the  output  matrix  C  is  represented  as 


Co 


1  0 


(7.22) 
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Now,  by  way  of  Equations  7.20,  7.21,  and  7.22,  the  HBSSM  presented  in  Equation 
5.43  can  now  be  created  for  the  linear  time  periodic  flapping  blade  model  of  the  form 
presented  in  Equation  7.23 

x  =  (A  —  JV)x  +  Bu 

y  =  Cx  +  Vu  (7.23) 

where  the  matrix  (A  —  Af),B,C,T>  of  the  HBSSM  are  in  block  Toeplitz  form. 

Now  that  a  single  flapping  blade  model  is  presented  in  HBSSM,  as  in  Equation 
7.23,  the  entire  rotor  system  can  now  be  modeled  by  combining  multiple  blade  models 
as  one  state  space  system.  As  noted  in  the  beginning  of  this  Section,  a  four  bladed 
rotor  system  is  desired.  To  begin  with,  four  individual  flapping  blades  represented 
each  in  the  form  of  equations  7.16,  7.18,  and  7.18  will  be  used  to  form  the  basic  rotor 
model,  as  seen  below. 


x(ijj)  =  A{^i)x{^})  +  B{ij})u{ij}) 

yW  =  C'OMVO  (7-24) 


where  the  plant  matrix,  A(i]S)  of  the  combined  system  is  represented  as 


[Ai\ 

[^2] 

[M 

[M 


(7.25) 
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and  the  control  matrix,  B(i]i)  of  the  combined  system  is  represented  as 


[Si] 

m 

[S3] 

[S4] 


(7.26) 


where  the  off  diagonal  terms  in  the  above  Equations  are  null.  As  a  note,  the  subma¬ 
trices  An,  n  =  1 ...  4  of  Equation  7.25  are  identical  to  each  other  and  are  individually 
equal  to  A (if;)  of  the  single  blade  Equation  7.18.  The  same  is  true  for  the  submatrices 
Bn,  n  =  1 . . .  4  of  Equation  7.26  which  are  individually  equal  to  B(ijj)  of  the  single 
blade  Equation  7.18.  Since  this  system  contains  no  blade  lag,  the  flapping  motion 
of  each  blade  is  uncoupled  from  any  other  blade.  As  this  system  now  contains  four 
blades,  each  with  states  (3,f3  and  having  control  input  0con,  the  state  and  control 
matrices  must  now  be  redescribed  as 


x(ijj) 
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As  a  note,  the  subscript  indices  in  Equations  7.27  and  7.28  refer  to  the  individual 
blade  with  which  they  are  associated. 

In  order  to  reconstruct  the  HBSSM  for  a  four  bladed  rotor,  the  Fourier  coeffi¬ 
cients  represented  in  Equations  7.20,  7.21,  and  7.22  will  be  expanded  to  contain  all 
four  blades  instead  of  a  single  blade.  For  ease  of  representation,  each  matrix  from 
the  HBSSM  that  is  associated  with  a  specific  blade  will  be  denoted  by  a  subscript, 
numbered  by  the  number  of  the  blade.  For  example  the  zero  harmonic  contribution 
of  the  second  blade  of  the  A  matrix  will  be  represented  as  Aq2.  As  all  of  the  blades 
are  assumed  identical  the  Fourier  series  coefficients  represented  in  Equations  7.20, 
7.21,  and  7.22  for  a  single  blade  are  simply  repeated  for  every  blade.  Thus,  for  the 
example  of  the  second  blade,  A02  =  A0  of  Equation  7.20.  This  procedure  is  identical 
for  the  B  and  C  Fourier  matrix  formulations.  Using  this  rationale,  the  Fourier  series 
coefficients  of  the  four  bladed  rotor  system  are  presented  in  Equations  7.29,  7.30,  and 
7.31  below.  First,  for  the  A  matrix  Fourier  coefficient  representation: 
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Next,  for  the  B  matrix  Fourier  coefficient  representation: 


Bo 


Bx 


B2 


[B0l]  0  0  0 

0  [BQ2]  0  0 

0  0  [B03]  0 

0  0  0  [Bo4] 

[Bu]  0  0  0 

0  [Bl2]  0  0 

0  0  [Bl3]  0 

0  0  0  [Bu] 

[B2l]  ooo 
0  [B22]  0  0 

0  0  [B2a]  0 

ooo  [B24\ 


B- 1  =  B* 


2  =  B* 


(7.30) 


Finally,  the  C  matrix  Fourier  coefficient  representation  contains  only  the  zero  har¬ 
monic  as  it  is  not  periodic  as  stated  above. 


[Oh]  ooo 
o  [Co,]  o  0 

o  0  [Co,]  0 

0  0  0  [CoA] 


(7.31) 


Now,  using  the  same  approach  for  the  HBSSM  as  for  the  single  blade,  the  HBSSM  of 
the  four  bladed  rotor  system  can  now  be  created.  The  Fourier  coefficient  matrices  for 
the  four  bladed  system  represented  in  Equations  7.29,  7.30,  and  7.31  are  used  to  form 
the  block  Toeplitz  matrices  (A  —  A f),  B ,  C,  V  used  to  create  the  four  bladed  HBSSM 
of  the  form  listed  in  Equation  7.23. 

With  the  completion  of  the  four  bladed  rotor  system  in  HBSSM  form,  we  are 
able  to  now  able  to  perform  Cramer-Rao  bound  analysis  on  the  individual  system  pa- 
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rameters  developed  from  system  identification.  Furthermore,  we  can  now  develop  an 
optimal  LQR  controller  to  attenuate  the  out  of  plane  rotor  vibrations.  The  structure 
of  the  controller  will  first  be  explained. 

7.f  Design  of  a  Optimal  Vibration  Attenuation  Controller  for  a  Linear 
Time  Periodic  Rotor  System 

The  intent  of  this  research  is  to  reduce  to  a  suitable  level  the  out  of  plane  vibra¬ 
tions  of  the  main  rotor  system.  However,  due  to  the  modeling  constraints  imposed 
by  existing  control  methods,  several  prerequisites  to  controller  development  must  first 
be  addressed  before  the  control  development  can  go  forward.  The  first  prerequisite 
is  namely  that  the  model  must  be  represented  as  LTP  for  accuracy.  The  second  pre¬ 
requisite  is  the  formation  of  a  state  space  representation  of  the  LTP  rotor.  These 
prerequisites  have  been  addressed  in  the  previous  chapters,  and  thus  an  adequate 
vibration  attenuation  controller  can  now  be  explored  and  developed. 

As  stated  in  chapter  III,  several  control  methodologies  currently  exist  in  the  area 
of  rotorcraft  control,  however  optimal  control  [7,36,40]  will  be  the  focus  of  this  research 
due  to  the  robust  characteristics  to  systems  with  poor  state  knowledge.  Furthermore, 
linear  optimal  methods  were  selected  for  use  as  they  inherently  use  parametric  models 
for  the  controller  development.  This  is  especially  attractive  as  the  one  of  the  key 
aspects  of  this  research  is  to  both  identify  and  verify  a  the  parameters  of  a  parametric 
linear  model.  Finally,  the  guaranteed  stability  margins  of  infinite  upward  and  one 
half  downward  gain  margins  are  necessary  for  a  system  that  is  subject  to  poor  system 
parameters.  The  next  Section  will  briefly  cover  the  development  of  a  LTP  linear 
quadratic  regulator  followed  by  a  tracking  regulator  design  used  to  eliminate  rotor 
vibrations. 

7.5  The  Linear  Time  Periodic  Linear  Quadratic  Regulator 

The  linear  quadratic  regulator  (LQR)  is  a  well  established  method  for  optimal 
control  for  a  linearized  plant  model.  The  LQR  method  is  built  upon  optimising  a 
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quadratic  cost  function  while  assuming  perfect  state  knowledge.  As  the  LQR  con¬ 
troller  has  been  adequately  described  for  both  LTI  [2]  and  LTP  [42]  systems,  a  brief 
description  of  the  regulator  will  be  presented.  The  formulation  of  the  linear  quadratic 
regulator  is  similar  for  both  LTI  and  LTP  systems  by  the  relation  of  the  harmonic 
balance  state  space  model.  As  such,  the  presentation  of  LQR  will  begin  in  terms  of 
the  LTI  system  and  then  present  the  LTP  dual. 

The  model  representing  the  linear  state  equations  for  the  LQR  are  in  state  space 
form,  as  presented  in  Equation  7.32. 

x(t)  =  Ax(t)  +  Buu(t )  (7.32) 

where  Bu  is  the  controllable  input  matrix.  The  quadratic  cost  function  J  that  forms 
the  core  of  the  LQR  is  presented  in  Equation  7.33. 

J(x(t),u(t))  =  \xT (t f)Hx(t f)  +  \  f*f  [xT(t)Q,x(t)  +  uT (t)Hu(t)~\  dt  (7.33) 

where  the  matrices  Q,H  are  positive  semidefinite  state  weighting  matrices  and  the 
matrix  R  is  a  positive  definite  control  weighting  matrix.  The  Riccati  equation  is 
required  to  provide  the  optimal  control  solution  ,  as  it  provides  the  minimum  to 
Equation  7.33.  A  brief  review  of  the  Riccati  equation  will  now  be  presented,  as  it 
provides  the  key  link  between  an  LTI  and  LTP  linear  quadratic  regulator. 

7.5.1  The  Riccati  Equation.  The  most  common  method  of  solving  the 
constrained  quadratic  cost  associated  with  the  LQR  method  is  the  use  of  the  Riccati 
equation.  The  Riccati  equation  is  a  nonlinear  matrix  differential  equation  which 
solves  for  the  matrix  of  proportionality  P(t)  between  the  constraint  costates  and 
system  states  by  direct  integration  backward  in  time.  This  is  possible  as  the  Riccati 
equation  has  only  final  conditions. 
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The  solution  to  the  optimal  control  problem  is  essentially  reduced  to  a  solving 
for  P(t),  as  the  optimal  control  law  is  defined  as 

u(t)  =  K(t)x(t ) 

u(t)  =  - R~lBlP{t)x{t )  (7.34) 

where  K(t)  is  the  time  varying  gain  matrix  for  LTI  systems.  Thus,  —  R^1B^P{t)  = 
K(t )  .  Furthermore,  the  matrices  Q  and  R  are  from  the  cost  function  of  Equation 
7.33.  Additionally,  the  state(s)  x(t)  are  assumed  to  be  from  perfect  and  total  state 
knowledge. 

The  value  of  P(t)  associated  with  the  optimal  solution  is  found  by  integrating 
the  Riccati  equation,  represented  in  Equation  7.35  backward  in  time  from  the  final 
condition,  P (tf). 

P{t)  =  —P(t)A  -  ATP(t)  -  Q  +  P(t)BuR~lBlP(t )  (7.35) 

where  the  final  condition  is  given  as 

P(t,)  =  H.  (7.36) 

The  optimal  control  law  is  formed  by  inserting  the  result  of  Equation  7.35  into  Equa¬ 
tion  7.34. 

The  formulation  of  the  Riccati  equation  in  Equation  7.35  is  valid  for  a  LTI 
system.  This  formulation  can  be  expanded  for  a  LTP  system  of  the  form  presented 
in  Equation  7.23  by  the  use  of  the  harmonic  balance  state  space  transformation. 
Thus,  using  the  LTP  HBSSM  state  model,  the  LTP  Riccati  equation  is  represented 
as  described  in  reference  [42] 

0  =  P 2  (*^2  ~~  A/^)  +  (A.2  —  0-2  —  V2B2uTZ2  lBlV2  (7.37) 

(7.38) 
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where  the  matrices  A2,  B2u,  V2,  Q2  and  1Z2  are  represented  as 


Re  (A)  Im(A) 


(7.39) 


The  formulation  of  the  remaining  matrices  B2,V2,  Q 2  and  1Z2  are  performed  in  an 
identical  manner  as  that  seen  in  Equation  7.39.  The  computation  of  J\f2  is  also  per¬ 
formed  in  a  similar  manner  as  that  of  Equation  7.39,  with  the  note  that  Re(A7)  =  0. 
Thus, 


0  -jA 
~jA  0 


A/2 


(7.40) 


The  augmented  matrices  for  the  LTP  Riccati  equation  are  necessary  to  transform  the 
equation  from  a  complex  algebraic  representation  to  a  real  algebraic  representation. 
This  transformation  makes  available  the  many  algebraic  Riccati  equation  solvers, 
which  is  advantageous  as  they  are  both  numerically  stable  and  accurate. 

The  solution  to  the  optimal  gain  in  the  LTP  case  retains  the  same  form  to  that 
of  the  LTI  case,  with  the  exception  that  the  LTP  matrices  are  in  Toeplitz  form.  . 
Thus,  by  inserting  the  transformed  Toeplitz  matrices  B2U,V2,  and  IZ2  into  Equation 
7.34,  the  LTP  gain  is  presented  as 


u{t)  =  —K(t)x{t) 
u(t)  = 


(7.41) 


where  /C(t)  is  the  time  periodic  gain  matrix.  As  with  the  Toeplitz  matrices  A,  B,C,  andT> 
of  the  HBSSM,  the  time  periodic  gain  matrix  A 2(t)  is  of  doubly  infinite  dimension. 

For  practical  application,  the  gain  matrix  must  be  truncated  in  order  to  execute  on  a 
digital  computer.  This  will  be  further  covered  in  Section  7.7.3. 
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7.6  Vibration  Reduction  via  Reference  Tracking  Control 

The  linear  quadratic  regulator,  which  is  designed  to  drive  the  states  of  the 
system  to  zero,  must  be  augmented  to  act  as  a  vibration  controller.  A  controller,  as 
compared  to  a  regulator,  can  track  a  reference  input  by  forcing  the  reference  input  to 
match  the  reference  output.  Thus,  the  LQR  can  be  modified  to  a  Reference  Tracking 
Controller  by  way  of  augmenting  the  LTP  system  model  with  an  additional  error  state, 
which  is  simply  the  error  between  the  reference  input  and  output.  This  additional 
state  is  handled  by  the  reference  tracking  method,  integral  feedback  control. 

Integral  feedback  is  a  method  in  which  an  error  state,  e  is  used  to  zero  out 
errors  between  a  constant  reference  input  signal  and  the  reference  output  from  the 
plant.  The  reference  signal,  in  the  case  of  this  controller,  will  be  a  blade  flapping 
angle,  /3ref,  and  will  be  covered  in  greater  detail  in  Section  7.6.1. 

The  error  state  is  simply  the  integral  of  the  differential  equation,  denoted  in 
Equation  7.42 


r  -  y(t) 
r  -  Cmx(t) 


e 


(7.42) 


where  r  is  the  reference  input  signal  and  Cm  is  the  measurement  output  matrix.  As 
a  note,  Cm  =  Cy,  where  Cy  is  the  state  output  matrix.  The  error  state  is  added  to 
the  rotor  system  model  to  form  the  augmented  system  seen  in  Equation  7.43. 


x(t)  =  Ax(t )  +  Buu(t)  +  Brr 


(7.43) 
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where  the  expanded  plant,  A,  control  Bu,  and  reference  Br  matrices  are  expanded  as 
seen  in  Equation  7.44 


x(t) 

— 

a  ; 

0 

x{t) 

+ 

Bu 

u(t)  + 

0 

e(t) 

—Cm  : 

0 

e(t) 

0 

I 

(7.44) 


The  control  law  u(t) 
state,  as  seen  below. 


—K{t)x  is  augmented  to  reflect 


the  addition  of  the  error 


u{t) 


e(t) 


Kx(t)  i  Ke(t) 


x(t) 

e(t) 


(7.45) 


The  tracking  controller  is  formed  by  inserting  Equation  7.45  into  7.44,  as  seen  in 
Equation  7.46. 


x{t) 

A 

0 

x(t) 

Bu 

x(t) 

= 

+ 

-Kx(t)\  —  Ke(t) 

.  . 

-crn 

0 

.  e W  . 

0 

.  . 

(7.46) 

The  above  representation  of  the  tracking  controller  is  further  simplified  by  multiplying 
out  the  gain  and  control  matrices,  as  seen  below. 


x(t) 

— 

A 

0 

+ 

- BuKx(t )  : 

— BuKe(t ) 

x{t) 

.  . 

-cm 

0 

0 

0 

.  e . 

(7.47) 
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Figure  7.3:  Reference  Tracking  Controller  Block  Diagram. 
The  tracking  controller  is  finalised  by  combining  like  terms 


x(t) 

A  -  BuKx(t )  : 

-BuKe(t) 

x(t) 

0 

= 

+ 

e(t)  _ 

—  C  : 

0 

.  e . 

I 

The  reference  tracking  controller  described  above  is  seen  in  Figure  7.3. 

The  reference  tracking  controller  described  so  far  is  applicable  to  a  LTI  system. 
This  is  easily  overcome  by  using  the  existing  LTP  optimal  gain  calculation  method 
covered  in  Section  7.5,  with  the  exception  of  using  an  augmented  system  with  an 
error  state.  Thus,  using  a  process  similar  to  that  covered  in  Equations  7.29  through 
7.31  the  Fourier  series  components  of  the  A,B,  and  C  matrices  can  be  assembled  with 
the  addition  of  an  error  state,  reference  tracking  controller  can  be  adapted  to  a  LTP 
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system.  Therefore,  the  augmented  state  vector  is  as  follows 


x(ijj) 


Pi 

Pi 

P2 

P2 

Ps 

Ps 

Pa 

Pa 


(7.49) 


which  is  the  same  as  Equation  7.27  with  the  exception  of  the  error  state.  As  a  note, 
the  control  input  vector  is  the  same  as  Equation  7.28.  Now,  for  the  augmented  matrix 
A,  the  Fourier  coefficient  representation  is  presented  with  the  additional  error  state 
as 
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Ao 


Ai 


A2 


[A0l]  0  0  0  0 

0  [i40a]  0  0  0 

0  0  [Aj3]  0  0 

0  0  0  [AoJ  0 

—  [Co]  _  [Co]  “[Co]  “[Co]  0 

[At^  0  0  0  0 

0  [Ala]  0  0  0 

0  0  [Al3]  0  0 

0  0  0  [Au]  0 

0  0  0  0  0 

[A2l]  0  0  0  0 

0  [A22]  0  0  0 

0  0  [A23\  0  0 

0  0  0  [A2i]  0 

0  0  0  0  0 


A~-!  =  A\ 


A- 2  — 


An 


(7.50) 
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Next,  for  the  augmented  matrix,  B  the  Fourier  coefficient  representation  is  as  follows 


B0 


Bo 


[Bq,]  0  0  0 

0  [Bq2]  0  0 

0  0  [B03]  0 

0  0  0  [B04] 

0  0  0  0 
[Bu]  0  0  0 

0  [Bl2]  0  0 

0  0  [Bl3]  0 

0  0  0  [Bu] 

0  0  0  0 
[b2i]  0  0  0 

0  [B22]  0  0 

0  0  [B23]  0 

0  0  0  [B2a] 

0  0  0  0 


B- !  =  B{ 


(7.51) 


As  before,  the  C  matrix  Fourier  coefficient  representation  contains  only  the  zero 
harmonic  as  it  is  not  periodic  as  stated  above.  It  does  contain,  however,  the  reference 
output  as  the  sum  of  all  four  blade  flapping  angles,  (3 .  The  augmented  matrix  C  is 
not  required  for  the  gain  calculation,  however,  it  will  be  needed  in  the  development  of 
the  reference  tracking  controller.  The  selection  of  the  reference  output  will  be  covered 
in  Section  7.6.1. 


Co 


[C0l]  0  0  0  0 

0  [C02]  0  0  0 

0  0  [C03]  0  0 

0  0  0  [Co,]  0 

c„,  [Co,]  [C„,]  [Co.]  0 


(7.52) 
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Now  that  the  Fourier  representations  of  the  state  and  control  matrices  have  been 
assembled  the  block  Toeplitz  forms  A  and  B  can  be  created  using  the  method  covered 
earlier.  From  this  point,  the  LQR  gains  can  be  computed  in  an  identical  fashion 
as  covered  in  Section  7.5  using  these  matrices.  The  remainder  of  the  LTP  tracking 
controller  development  is  identical  to  that  of  an  LT1  system. 

7.6.1  Tracking  Control  Applied  to  Out  of  Plane  Rotor  Vibrations.  In  order 
to  use  the  tracking  controller  developed  in  Section  7.4  an  adequate  reference  input 
and  output  must  be  defined.  As  noted  in  chapter  II,  the  out  of  plane  vibrations  in  a 
helicopter  rotor  system  are  caused  by  differences  in  blade  lift  when  each  blade  in  the 
rotor  system  is  considered.  In  this  work,  each  blade  is  considered  to  have  identical 
profile  and  thus  if  each  blade  has  an  identical  flap  track  over  one  period,  the  lift  of  each 
blade  is  then  identical.  Thus,  an  input  reference  signal  is  generated  by  considering 
the  combined  blade  flapping  angles  of  a  perfect  rotor.  This  is  represented  in  Figures 
7.4  and  7.5,  which  is  representative  of  the  blade  flapping  induced  by  an  input  of  one 
degree  pitch  0con  i_4  =  1  for  each  blade  at  an  inflow  ratio  of  fi  —  .3.  The  individual 
blade  flapping  angles  of  all  blades  in  the  rotor  system  are  summed  at  each  time  step, 
as  presented  in  Figure  7.5.  This  value  is  the  optimal  reference  value,  as  in  this  case 
all  the  blade  lifts  are  identical,  thus  having  no  out  of  plane  vibrations.  The  output 
reference  value  is  the  summation  of  each  of  the  blade  flapping  angles,  as  seen  in  the 
output  matrix  C  in  Equation  7.52.  The  steady  state  value  presented  in  Figure  7.5  is 
the  reference  input  that  will  be  applied  to  the  test  rotor  system  . 

The  performance  of  the  vibration  controller  is  evaluated  in  Section  7.8  by  in¬ 
ducing  a  lift  imbalance  in  the  rotor  system  an  then  reviewing  the  ability  of  the  system 
to  eliminate  the  ensuing  vibration.  To  induce  the  lift  imbalance,  pitch  inputs  of  one 
or  more  of  the  blade  pitch, @con,  differing  from  a  the  desired  setting  of  one  degree 
will  be  set.  This  will  generate  an  imbalance  in  lifts  between  the  individual  blades, 
which  ultimately  will  cause  an  out  of  plane  oscillatory  vibration.  The  reference  track¬ 
ing  function  of  the  vibration  controller  will  adjust  the  blade  pitches  of  each  blade 
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Blade  Flapping  Angle, 


Blade  Flapping  Angle  For  Blades  1-4 


Figure  7.4:  Individual  Flap  Angles  of  the  Four  Blades. 
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Combined  Blade  Flapping  Angle, 


Difference  in  Blade  Flapping  Angles 


Figure  7.5:  Combined  Flap  Angles  of  the  Four  Blades. 
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until  the  input  and  output  references  are  the  same,  thus  eliminating  the  rotor  vibra¬ 
tions.  More  on  this  will  be  covered  in  Section  7.8,  which  details  an  application  of  this 
controller  on  a  rotor  system  that  has  incorrect  pitch  input. 

The  next  Section  will  outline  the  application  of  this  controller,  while  demon¬ 
strating  how  the  Cramer-Rao  bound  can  validate  the  control  solution. 

7. 7  Controller  Model  Parameter  Validation  via  the  Cramer-Rao  Bound 

An  accurate  vibration  controller  has  a  fundamental  underpinning;  an  accurate 
model  which  the  controller  is  based  upon.  This  is  important,  as  the  controller  will  only 
be  as  effective  as  the  model  is  accurate.  For  an  adaptive  controller,  such  as  one  that 
identifies  the  controller’s  model  ’on  the  fly’,  a  method  such  as  the  Cramer-Rao  lower 
bound  that  can  validate  the  model  will  in  turn  ensure  accurate  control.  This  Section 
will  develop  the  Cramer-Rao  lower  bound  to  validate  model  parameters  which  are 
based  upon  system  identification.  The  Cramer-Rao  bounds  will  then,  in  Section  7.8, 
be  used  to  validate  the  effectiveness  of  an  adaptive  version  of  the  vibration  controller 
described  in  Section  7.6. 

This  Section  will  begin  by  describing  the  calculation  of  the  Cramer-Rao  bound 
for  parameters  of  the  rotor  system  model.  The  Section  will  conclude  by  showing  the 
Cramer-Rao  bounds  for  a  select  model  parameter,  with  the  intent  to  demonstrate  the 
validity  of  model  parameters  developed  from  an  online  system  identification  method. 

7.7.1  Cramer-Rao  Bound  Calculation  for  a  LTP  Rotor  Model.  As  stated 
in  Section  7.7,  the  accuracy  of  the  individual  parameters  that  make  up  a  model,  in 
this  case  a  HBSSM  rotor  model,  is  essential  to  developing  an  accurate  vibration  con¬ 
troller.  In  the  case  of  poorly  identified  parameters  from  system  identification,  the 
model  developed  from  those  parameters  will  lead  to  inaccurate  LQR  gains.  This 
will  ultimately  lead  to  poor  controller  performance  in  terms  of  tracking  error  or  ex¬ 
cessive  control  input  to  compensate  for  modeling  errors.  This  Section  will  develop 
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the  Cramer-Rao  bounds  needed  to  detect  poorly  identified  system  parameters  in  an 
attempt  to  provide  better  online  controller  performance. 

The  development  of  the  Cramer-Rao  lower  bound  for  the  LTP  rotor  model 
described  in  Section  7.3  will  begin  by  recalling  the  work  covered  in  Section  6.6.2.  For 
starters,  the  Cramer-Rao  lower  bound  can  be  computed  by  taking  the  square  root 
of  inverse  of  the  block  Toeplitz  Hessian,  7d,  as  seen  below  in  Equation  7.53.  As  a 
note,  the  index  j  relates  to  the  parameter  for  which  the  Cramer-Rao  bound  is  being 
developed. 


CK,~  JvrXj  (7.53) 

The  Hessian  is  formed  by  summing  the  square  of  the  partial  derivative  of  the  frequency 
response  function,  ^ ,  with  the  measurement  signal  noise  weighting  matrix,  W, 
over  all  input  frequencies,  ujf.  The  partial  derivative  frequency  response  function 
is  described  using  the  state  space  matrices  A,  A f,B,  C,  and  V  from  Eq.  7.23  in  a 
modified  form  of  Equation  7.55,  as  seen  below. 

|  ■  AnNylB  +  ®i 

BA 

=  C(tUfI  -  AHU)~1  —  {LwfI  -  AHNylB 

c)B 

+  C{aufI  -  AHN)-1  — 

where,  as  before,  the  matrix  AHM  is  defined  as 

AHAf  =  A  -  Af  .  (7.56) 

The  above  Equations  7.55  and  7.56  redefine  the  partial  derivative  of  the  response  error 
from  Equation  6.27  to  an  expression  containing  a  doubly  infinite  block  of  complex 
Fourier  coefficients  in  Toeplitz  form.  To  complete  the  computation  of  the  Cramer- 
Rao  lower  bound  for  LTP  systems  defined  in  block  Toeplitz  form,  Equation  6.48  is 


(7.54) 

(7.55) 
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used  to  in  Equation  7.57  to  define  the  Hessian,  77,  in  a  manner  analogous  to  Equation 
6.31. 


1=1 


dQ j,  k{^  fi) 

dp 


(7.57) 


In  Equations  7.20  through  7.22  the  Fourier  series  coefficients  matrices  of  an 
individual  rotor  blade,  Equation  7.18  and  7.18,  were  presented.  As  covered  previously, 
these  matrices  are  used  to  form  the  block  Toeplitz  matrices  A  and  B.  The  partial 
derivatives  ^  and  ||  are  developed  for  Equation  7.55  in  an  identical  manner  by  using 
the  definitions  for  the  Fourier  coefficients  '-A  and  as  seen  in  Equations  7.58  and 
7.59.  It  is  important  to  note  that  in  this  case,  the  parameter  p  with  which  the  partial 
derivative  is  taken  in  the  Lock  number ,  7.  The  rationale  for  this  parameter  selection 
will  be  covered  in  Section  7.8. 
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As  before,  by  developing  the  Fourier  series  coefficients,  the  block  Toeplitz  forms  of 
the  plant  and  control  matrices  were  then  populated.  Similarly,  the  same  method  is 
applied  to  form  the  partial  derivatives  ^  and  ^  by  using  Equations  7.58  and  7.59 
to  form  Equation  7.60. 
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(7.60) 


where  is  performed  in  a  likewise  manner. 

The  necessary  components  of  Equation  7.53  are  now  complete  and  the  Cramer- 
Rao  bound  can  be  calculated  for  the  identified  Lock  parameter,  7,  of  the  blade  model. 
The  next  Section  will  show  results  of  the  Cramer-Rao  bound  for  the  parameter  7, 
which  is  derived  from  the  frequency  domain  system  identification  method  developed 
by  ffwang  [12]  for  LTP  systems. 


1.1.2  Cramer-Rao  Bounds  of  Lock  Number,  7.  In  this  Section  the  Lock 
number,  7,  of  the  rotor  blade  model  represented  by  Equations  7.18  and  7.18  will 
be  validated  for  accuracy  by  way  of  the  Cramer-Rao  bound.  The  parameters  to  be 
validated  are  first  collected  by  way  of  a  frequency  domain  system  identification  method 
developed  by  Hwang  [12]  for  LTP  systems  over  a  range  of  input  frequency,  a and 
measurement  noise,  Sv.  Once  these  values  are  obtained  for  all  frequency  inputs  and 
noise  variances,  the  corresponding  Cramer-Rao  bound  calculated  by  way  of  Equation 
7.53  is  then  compared  to  the  identified  value. 
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The  system  identification  process  begins  by  using  a  blade  model  consisting  of 


the  following  parameters,  as  defined  in  table  7.61 


Advance  Ratio,  /i 

Lock  Number,  7 

H  =  0.3 

7  =  8 

The  system  identification  process  proceeds  by  inputing  an  oscillation  of  frequency 
Uf  in  ©con-  Each  system  identification  run  is  performed  at  a  single  input  frequency, 
where  the  range  of  Uf  =  u)p(0,  0.05.  0.1,  0.2,  0.3,  0.4,  0.5),  where  up  =  1.  As  a 
note,  both  the  system  identification  method  and  the  Cramer-Rao  bound  method  can 
accept  a  range  of  multiple  input  frequencies  for  each  run,  such  as  a  signal  chirp.  This 
method  was  not  selected  for  this  work,  however,  as  the  parameter  validity  was  desired 
at  individual  frequencies  so  that  a  range  of  acceptable  parameters  could  be  identified. 
This  will  be  addressed  as  the  plots  of  the  Cramer-Rao  bounds  are  discussed. 

As  each  system  is  stimulated  at  one  of  the  values  of  u>f,  the  system  output 
measurement  is  corrupted  by  white  noise  having  spectral  density,  Sv.  In  this  work, 
a  range  of  spectral  densities  Sv  =  1,  2,  3,  4  was  considered.  For  all  values  in  the 
range  of  Sv,  a  run  is  made  for  every  input  value  in  the  range  of  07.  Figures  A.l,  A. 2, 
A. 3,  and  A. 4  of  Appendix  A  are  plots  of  the  identified  value  of  7  at  each  ujf  with  the 
corresponding  Cramer-Rao  bound  overlaid.  Each  plot  is  representative  of  the  one  of 
the  values  of  Sv  in  the  defined  range  above. 

By  reviewing  the  comparison  plots  mentioned  above  in  Appendix  A,  it  is  clear 
that  both  input  frequency,  07,  and  the  intensity  of  the  measurement  noise,  Sv,  directly 
affect  the  quality  of  the  estimate  of  the  parameter,  7.  First,  consider  the  effect  of  the 
input  frequency,  a 7,  on  the  accuracy  of  the  parameter  estimate.  As  was  discussed  in 
chapter  VI,  the  Cramer-Rao  bound  presents  a  scaled  inverse  of  the  frequency  response 
function.  Thus,  as  frequency  response  falls  off  the  Cramer-Rao  bound  begins  to  grow 
in  magnitude.  This  is  seen  in  Figure  7.6,  which  plots  the  Cramer-Rao  bound  for  the 
blade  model  for  all  input  frequencies  from  0  <  07  <  0.507.  As  an  example,  this 
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COMBINED  CRAMER-RAO  PLOT,  NOISE  VARIANCE  =  1-4 


Figure  7.6:  Cramer-Rao  Bounds  for  all  Frequencies,  Uf,  all  Noise  Values,  Sv. 

plot  reveals  that  at  low  input  frequencies  the  corresponding  Cramer-Rao  bound  has 
a  large  magnitude.  This  is  due  to  the  poor  frequency  response  of  the  blade  model 
at  low  frequency  input.  Secondly,  this  plot  also  presents  the  Cramer-Rao  bound 
over  this  range  for  each  value  of  Sv,  which  demonstrates  the  effect  of  measurement 
signal  noise  on  the  magnitude  of  the  bound.  When  considering  the  computation  on 
the  Hessian,  as  seen  in  Equation  7.57,  a  weighting  factor  W  scales  the  effect  of  the 
partial  derivative  of  the  frequency  response  function.  As  a  reminder,  this  weighting 
factor  is  formed  by  taking  the  inverse  of  the  prediction  error  covariance  matrix,  R.  As 
an  additional  note,  the  prediction  error  covariance  matrix  in  this  case  is  simply  equal 
to  the  power  spectral  density  of  the  measurement  noise,  Sv.  Tims,  W  =  i?”1  and 
R  =  Sv.  With  this  in  mind,  as  signal  noise  increases  the  magnitude  of  the  Hessian 
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decreases.  This  ultimately  increases  the  magnitude  of  the  Cramer-Rao  bound,  as  seen 
in  Equation  7.53.  This  is  expected,  as  any  degradation  in  measurement  signal  quality 
will  effect  the  input-output  relation  of  the  system.  The  varying  levels  of  noise  are 
depicted  in  the  plot  as  Var±,  Var2,  V ar3,  and  Var 4  depict  the  input  noise  levels, 
Sv  —  1,  Sv  —  2,  Sv  —  3,  and  Sv  =  4,  respectively. 

The  plots  of  Appendix  A  depict  the  relationship  between  the  identified  value 
of  7  at  each  Uf  and  Sv  with  the  corresponding  Cramer-Rao  bound  .  To  demonstrate 
this  relationship,  two  cases  will  be  discussed.  The  first,  as  seen  in  Figure  7.7  depicts 
Cramer-Rao  bounds  for  low  frequency  response  and  high  signal  noise.  The  second 
case,  as  seen  in  Figure  7.8  depicts  Cramer-Rao  bounds  for  high  frequency  response 
and  low  signal  noise.  The  analysis  will  begin  with  case  one. 

Case  one  considers  cases  of  poor  measurement  signal  quality  and  low  signal 
input  frequency,  such  as  the  case  of  Figure  7.7  at  a 7  =  0.  Here,  the  accuracy  of  7 
is  low,  as  the  identified  value  7  =  5.2,  where  as  the  true  value  of  7  =  8.  This  is 
expected,  however,  as  the  magnitude  of  the  Cramer-Rao  bound  is  almost  equal  to  9. 
As  the  Cramer-Rao  bound  is  simply  the  standard  deviation  of  the  expected  value  of 
7,  the  bound  indicates  the  expected  value  of  7  to  be  in  the  range  of  0.5  <  7  <  8.9. 
With  a  standard  deviation  this  large,  it  is  obvious  that  this  identified  value  of  7  is  of 
low  value  and  needs  to  be  re-evaluated  or  more  importantly,  discarded  altogether. 

The  opposite  is  true  for  case  two,  which  considers  values  of  7  in  Figure  7.8, 
evaluated  at  a if  =  0.5.  In  this  case,  both  the  measurement  noise  is  low  and  the 
input  frequency  is  high,  thus  corresponding  to  a  small  Cramer-Rao  bound.  As  would 
be  expected,  the  identified  value  of  7  is  nearly  perfect,  with  7  =  7.99.  This  cor¬ 
responds  to  a  standard  deviation  of  the  expected  value  of  7  to  be  in  the  range  of 
7.9  <  7  <  8.1.  In  this  case,  it  is  advisable  to  accept  this  identified  parameter  as 
accurate  as  the  bounds  are  small. 

As  discussed  in  Chapter  HI,  the  Cramer-Rao  bound  is  superior  to  traditional 
data  scatter  analysis  to  evaluate  the  validity  of  the  an  estimated  parameter.  This  was 
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Figure  7.8:  Cramer-Rao  Bounds  for  all  Frequencies,  Uf,  at  Noise  Value,  Sv  =  1. 
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demonstrated  in  Figure  3.2,  which  showed  that  tight  data  scatter,  a  common  used 
method  to  qualify  data  as  accurate,  does  not  directly  relate  to  an  accurately  estimated 
parameter.  Only  by  reviewing  the  magnitude  of  the  Cramer-Rao  bound  versus  the 
data  scatter  can  the  parameter  estimates  truly  be  determined  as  accurate.  There¬ 
fore,  to  demonstrate  the  usefulness  of  the  Cramer-Rao  bound  to  validate  parameter 
estimates,  the  Cramer-Rao  bounds  represented  in  Appendix  A  were  superimposed  on 
estimates  of  the  blade  Lock  number,  7.  The  data  scatter  represents  the  estimates  of 
the  Lock  number  from  100  individual  tests.  A11  example  of  this  analysis  is  presented 
in  Figures  7.9  and  7.10.  These  plots  present  the  Cramer-Rao  bounds  calculated  at 
each  the  input  frequencies,  a at  two  noise  levels,  Sv  =  1  and  Sv  =  4.  Upon  review 
of  both  Figures  it  is  clear  that  data  scatter  does  not  directly  relate  to  the  validity 
of  parameter  estimates.  Take  Figure  7.9  for  example.  By  review  of  the  values  of  7 
estimated  at  the  optimal  frequency  response  point,  ivj  =  0.5,  both  the  Cramer-Rao 
bound  and  the  data  scatter  agree.  However,  as  the  frequency  response  falls  off  as 
Uf  — *  0  the  Cramer-Rao  bound  and  the  corresponding  data  scatter  do  not  agree. 
This  is  most  pronounced  at  Uf  =  0.0,  where  the  data  scatter  is  very  tight,  however 
the  Cramer-Rao  bound  is  large.  If  one  was  to  evaluate  the  data  scatter  alone  at  this 
point  the  tight  data  scatter  may  lead  to  a  false  determination  that  the  parameter 
estimate  was  indeed  accurate.  Figure  7.10  can  be  evaluated  as  previously  described, 
however  this  plot  reveals  the  effect  of  greater  signal  noise  on  the  parameter  estimates. 
While  the  data  scatter  is  overall  less  tight  than  that  of  the  lower  noise  case  of  Figure 
7.9,  the  same  conclusion  can  be  made  regarding  the  false  readings  data  scatter  can 
produce  when  determining  parameter  estimates. 

The  size  of  the  bound  that  corresponds  to  a  poor  estimate  is  ultimately  left 
to  engineering  judgement,  as  to  what  is  acceptable  in  the  particular  application.  In 
the  next  Section  the  size  of  the  bound  will  be  evaluated  in  terms  of  the  accuracy 
of  the  vibration  suppression  controller  developed  in  this  chapter.  Before  this  is  cov¬ 
ered,  however,  a  brief  word  on  the  adaptation  of  the  HBSSM  for  the  use  in  a  digital 
computer  will  be  made. 
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Figure  7.9:  Cramer-Rao  Bounds  vs.  Data  Scatter  of  100  runs  for  all  Frequencies, 
Uf,  at  Noise  Value,  Sv  =  1. 
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IDENTIFIED  VALUES  OF  LOCK  PARAMETERy,  AND  COORESPONDING  CRAMER-RAO  BOUND 
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Figure  7.10:  Cramer- Rao  Bounds  vs.  Data  Scatter  of  100  runs  for  all  Frequencies, 

ujf,  at  Noise  Value,  Sv  =  4. 
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7.7.3  Determination  of  Model  Dimension.  The  Toeplitz  block  form  of 
Fourier  coefficients  is  key  to  developing  a  state  space  form  of  the  system  matrices 
that  define  a  model.  These  matrices,  as  described  previously,  are  considered  doubly 
infinite  as  they  contain  a  doubly  infinite  representation  of  the  Fourier  series.  This, 
however  is  untenable  for  application  in  code  development  on  a  digital  computer.  For 
code  development  to  take  place,  the  representative  Fourier  series  of  the  Toeplitz  trans¬ 
form  must  be  truncated  to  contain  N  harmonics.  To  determine  the  value  of  N,  the 
number  of  harmonics  was  increased  in  the  system  as  to  identify  convergence  in  the 
compared  system  responses.  Based  on  the  experience  of  Hwang  [12],  N=5  was  de¬ 
termined  to  produce  adequate  convergence  in  system  output  response.  The  same 
procedure  was  conducted  when  determining  the  number  of  harmonics  to  include  in 
the  LQR  gains,  as  they  are  also  Fourier  series  dependent.  The  closed  loop  system 
response,  as  will  be  seen  in  Section  7.8,  was  evaluated  for  convergence  while  increasing 
the  number  of  harmonics  in  the  LQR  gains.  In  this  case,  N=2  proved  adequate  for 
an  optimal  control  solution  based  on  response  convergence. 

7.8  Vibration  Controller  Validation  via  the  Cramer-Rao  Bound 

In  this  final  Section  of  the  chapter  the  performance  of  the  vibration  controller 
developed  in  Section  7.6  will  be  evaluated  based  on  the  accuracy  of  the  identified 
system  model  parameters,  in  this  case  7.  The  model  of  the  rotor  system  will  be  that 
of  Section  7.3,  with  the  distinction  of  having  three  blades  with  an  identical  blade  pitch 
input  QCon  =  1  deg,  while  one  blade  will  differ  with  0con  =  1.3  deg.  This  difference  in 
blade  pitch  input  will  generate  asymmetric  lift  in  the  rotor  system,  thus  causing  an  out 
of  plane  vibration  synonymous  with  that  of  the  ’track  and  balance’,  as  mentioned  in 
Section  7.6.1.  The  intent  of  this  Section  is  to  demonstrate  the  effect  poorly  identified 
system  parameters,  determined  by  the  Cramer-Rao  bound,  have  on  the  performance 
of  a  controller. 

The  effect  of  a  parameter  variation,  of  the  type  discussed  in  Section  7.8,  is  that 
the  model  in  which  the  parameter  resides  will  differ  from  the  true  system  it  is  to 
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represent.  This  presents  a  problem  when  this  model  is  used  in  controller  synthesis. 
In  this  case  the  optimal  control  solution  is  designed  based  upon  an  a  model  that  has 
error  due  to  incorrectly  estimated  parameters  when  compared  to  the  true  system. 
For  example,  Burl  [2]  states  that  most  control  system  designers  have  synthesised  a 
controller  with  good  performance,  simulated  the  controller  to  verify  the  performance, 
then  implemented  the  controller,  only  to  discover  that  the  performance  is  totally 
unacceptable.  These  discrepancies  between  the  mathematical  model  of  the  plant  and 
the  actual  system  are  denoted  as  model  perturbations.  In  this  case,  and  moreover  in 
general,  Burl  [2]  notes  that  uncertainty  in  the  plant  is  modeled  by  a  set  of  feedback 
perturbations.  In  these  cases,  a  robust  control  scheme  must  be  used  to  ensure  stable 
operation  to  overcome  the  modeling  limitations. 

As  was  described  earlier,  the  LQR  control  method  was  selected  as  the  controller 
synthesis  method  upon  which  the  vibration  reduction  system  is  based  upon.  This 
selection  was  made  based  on  the  relevance  to  the  problem;  the  rotor  model  used 
for  controller  synthesis  will  contain  bounded  uncertainty  in  the  plant  parameters,  in 
this  case  the  Lock  number,  7.  The  Lock  number  was  selected  as  it  is  a  key  factor 
in  determining  the  flap  response  to  the  blade  model  used  in  this  study.  The  linear 
quadratic  regulator  is  tolerant  to  feedback  perturbations  and  as  such  was  selected 
based  on  these  characteristics.  This  tolerance  is  due  to  the  guaranteed  bound  on  the 
smallest  destabilising  feedback  perturbation,  which  by  using  the  triangle  inequality 
[2,42]  provides  an  guaranteed  gain  margins  of  GM+  =  oo,GM~  <  |.  As  this  system 
contains  multiple  inputs,  phase  margins  are  irrelevant.  With  these  stability  margins 
the  LQR  is  recognised  as  having  robust  performance  for  feedback  perturbations,  as 
this  system  is  defined  as  having. 

The  ability  of  this  type  of  controller  to  deliver  adequate  reference  tracking  per¬ 
formance  while  having  feedback  uncertainty  may  not  necessarily  guarantee  overall 
adequate  system  performance.  This  is  because  the  controller  will  generally  require 
greater  control  input  to  compensate  for  the  modeling  discrepancies  mentioned  above. 
While  these  control  inputs  will  not  destabilise  the  system,  they  may  hinder  the  ability 
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of  the  system  to  perform  other  functions.  An  example  of  this  is  control  saturation, 
which  would  prohibit  the  pilot  of  a  helicopter  from  performing  other  controlled  ma¬ 
neuvers  if  the  vibration  controller  is  dominating  the  available  control  input  to  the 
system.  Another  example  would  be  that  excessive  control  input  may  drive  the  system 
to  points  of  aerodynamic  stall,  as  in  the  case  of  high  ©con  demands  by  the  control 
system.  The  following  Section  will  present  an  example  of  this  by  reviewing  the  per¬ 
formance  of  the  vibration  controller  subject  to  varying  cases  of  model  uncertainty,  in 
the  case  of  the  varying  parameter  7. 

7.8.1  Vibration  Controller  Performance  Evaluated  By  Cramer- Rao  Bound. 
This  Section  will  present  the  performance  of  the  vibration  controller  developed  in  this 
work  subject  to  parameter  inaccuracies.  The  system  will  be  reviewed  over  the  range 
of  7  as  defined  in  Section  7.7  where  the  Cramer- Rao  bound  of  the  parameter  was 
compared  to  the  corresponding  system  identified  parameter  value.  This  range  of  test 
conditions  will  allow  for  a  direct  comparison  of  the  controller  performance  in  terms 
of  control  demands  exhibited  by  the  controller  to  the  validity  of  the  parameter  7. 
This  comparison  will  reveal  that  while  the  vibration  controller  has  adequate  tracking 
performance,  the  control  requirements  may  be  unacceptable  by  the  terms  listed  above. 

The  vibration  controller  test  condition  is  defined  by  the  inflow  ratio  of  fi  =  .3 
and  Lock  number  7  =  8  for  the  actual  rotor  system  at  a  pumping  frequency  of 

up  =  1.  The  Q  and  R  matrices  which  define  the  LQR  controller  were  selected  as 
such;  Q  =  diag([0  0  0  0  0  0  0  0  2.0]);  where  the  weighting  on  each  blade  state  is  zero 
and  increased  weighting  on  the  error  state,  and  R  weights  equally  the  four  control 
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inputs  by  ^  on  the  diagonal, as  seen  in  Equations  7.62  and  7.63. 


Q 


0  0  0 
0  0  0 
0  0  0 
0  0  0 
0  0  0 
0  0  0 
0  0  0 
0  0  0 
0  0  0 


0  0  0 
0  0  0 
0  0  0 
0  0  0 
0  0  0 
0  0  0 
0  0  0 
0  0  0 
0  0  0 


0  0  0 
0  0  0 
0  0  0 
0  0  0 
0  0  0 
0  0  0 
0  0  0 
0  0  0 
0  0  2 


(7.62) 


R  = 


—  0 

100  w 

0 


0 

0 


0 

0 


1 

100 


0  — 
100 


0 

0 

0 


1 

100 


(7.63) 


0  0 

These  selections  of  weighting  matrices  were  made  to  demonstrate  good  vibration 
reduction  while  not  limiting  controller  input.  The  rotor  system  will  produce  an  out  of 
plane  vibration  by  defining  the  input  to  the  third  blade  to  have  a  bias  input  of  0.3  deg 
to  any  commanded  input.  The  remaining  blades  will  have  no  bias  and  will  produce 
the  desired  command  input.  This  bias  is  designed  to  replicate  an  incorrectly  adjusted 
pitch  linkage  on  blade  three  as  compared  to  the  remaining  blades  of  the  rotor  system. 
The  open  loop  blade  flapping  response  of  the  rotor  system  subject  to  the  defined  test 
conditions  is  presented  in  Figure  7.11.  This  unbalance  in  flap  angles  corresponds  to 
a  oscillation  in  the  summation  of  the  blade  flapping  angles  as  seen  in  Figure  7.12. 
The  oscillation  in  the  summation  of  the  blade  flapping  angles  seen  in  Figure  7.12  is 
indicative  of  a  out  of  plane  rotor  vibration  caused  by  asymmetric  lift  in  the  rotor 
system,  as  discussed  in  Section  7.6.1.  The  intent  of  vibration  controller,  as  outlined 
in  Section  7.6.1  is  to  smooth  the  blade  individual  blade  flapping  angles  such  that  the 
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Blade  Flapping  Angle, 


Blade  Flapping  Angle  For  Blades  1-4 


Figure  7.11: 


Blade  Flapping  Angles  For  an  Unbalanced  Rotor. 
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Combined  Blade  Flapping  Angle, 


Difference  in  Blade  Flapping  Angles 


Figure  7.12: 


Summation  of  the  Blade  Flapping  Angles  For  an  Unbalanced  Rotor. 
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oscillations  present  in  Figure  7.12  will  reduce  to  match  the  reference  input  signal  of 
Figure  7.4.  This  reference  signal  is  designed  to  reproduce  the  flapping  response  of  the 
four  identically  bladed  rotor  system,  having  parameters  7  =  8,^  =  0.3,  with  identical 
pitch  inputs,  ©con  =  1  deg,  as  described  in  Section  7.6.1. 

For  this  review,  the  test  controller  will  be  reviewed  as  28  individual  test  points, 
each  corresponding  to  a  specific  value  of  the  parameter  7  identified  at  a  specific  value 
of  input  frequency  a 7  and  measurement  noise  Sv.  As  a  note,  these  are  the  same 
parameters  as  were  reviewed  in  Section  7.7.2  and  thus  the  analysis  of  that  Section 
will  carry  over  here.  The  28  total  test  points  that  correspond  to  the  identified  values  of 
7  are  defined  by  the  input  frequencies,  07  =  o;p(0,  0.05.0.1,  0.2,  0.3,  0.4,  0.5),  and 
measurement  noise  spectral  densities  Sv  =  1,  2,  3,  4.  Once  again,  these  were  the 
same  parameters  used  to  identify  the  Lock  number  in  Section  7.7.2.  The  performance 
of  each  test  point  are  presented  in  Appendices  C,  D,  E,  and  F,  each  corresponding  to 
a  value  of  Sv  =  1,  2,  3,  4,  respectively.  A  review  of  the  plots  of  the  individual  flap 
angles,  calculated  LQR  gains,  control  usage,  and  vibration  controller  performance  in 
terms  of  tracking  performance  for  each  of  the  28  individual  test  cases  is  presented  in 
these  appendices. 

Of  these  28  cases,  two  cases  representing  the  best  and  worst  values  of  uij  and  Sv 
used  to  identify  7  will  be  compared  with  respect  to  the  Cramer- Rao  bound  to  establish 
a  connection  between  bound  magnitude  and  controller  effectiveness.  The  first  test  case 
is  representative  of  the  best  identified  value  of  7,  having  a 7  =  .5u;p  and  Sv  —  1.  These 
points  are  based  on  analysis  of  Section  7.7.2  when  correlating  the  Cramer- Rao  bound 
to  parameter  estimate  quality.  The  second  test  case  is  representative  of  the  identified 
parameter  7  having  lowest  quality,  which  occurs  at  Uf  =  0.0cnp  and  Sv  =  4. 

A  quick  review  of  the  tracking  performance  of  the  vibration  controllers  in  Figures 
7.13  and  7.14  for  the  best  and  worst  cases,  respectively,  reveals  little  difference  in 
overall  tracking  performance  of  the  reference  signal.  This  is  to  be  expected  based  on 
the  previous  discussions  on  the  robustness  of  the  LQR  to  perturbations.  This  does 
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Controller  Tracking  of  Reference  Input 


Controller  Tracking  of  Reference  Input,  Zoom  In 


Figure  7.13:  Best  Case:  Tracking  Performance  of  Vibration  Controller  for  case 

Uf  =  0.5u+,  Sv  =  1. 

not  indicate,  however,  that  both  control  designs  are  equivalent,  as  will  be  discussed 
next. 

A  further  review  of  the  input  control  required  to  achieve  this  level  of  vibration 
reduction  in  both  test  cases  is  more  revealing  in  term  of  the  differences  between  the 
two  cases.  The  required  inputs,  0con  for  blades  1-4  of  the  best  case  range  from 
approximately  —6.8  <  0con  <  +2.2  are  greatly  reduced  compared  to  the  worst  case, 
whose  corresponding  values  of  range  from  approximately  —9.3  <  0con  <  +3.3.  These 
results  are  seen  in  Figures  7.15  and  7.16,  which  once  again  represent  the  best  and  worst 
cases,  respectively.  The  difference  in  the  control  requirements  are  not  surprising,  as 
the  increased  model  perturbation  due  to  poor  identification  of  7  in  the  case  of  the  worst 
case  require  greater  control  input  to  compensate  for  what  is  essentially  an  unknown 
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Controller  Tracking  of  Reference  Input 


Controller  Tracking  of  Reference  Input,  Zoom  In 


Figure  7.14:  Worst  Case:  Tracking  Performance  of  Vibration  Controller  for  case 

tOf  =  O.OtUp,  Sv  =  4. 
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Figure  7.15:  Best  Case:  Control  Usage  for  case  Uf  =  0.5 u>p,  Sv  =  1. 
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Figure  7.16:  Worst  Case:  Control  Usage  for  case  u>f  =  0.0cup,  Sv  =  4. 
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model  disturbance,  as  discussed  previously.  This  disturbance  is  not  destabilising, 
however  it  has  a  tremendous  impact  on  the  availability  of  control  input.  This  test 
case  is  not  representative  of  a  real  system  ,  however  if  one  was  to  assume  that  the 
available  control  inputs  were  between  a  range  of  —10  <  @cori  <  10,  the  available 
control  authority  of  the  worst  case  would  be  almost  zero  in  terms  of  negative  ©con.  The 
best  system,  by  comparison,  still  has  adequate  control  authority  available  by  the  same 
standards.  Thus,  by  relating  the  Cramer-Rao  bound  to  the  parameter  perturbations 
in  a  manner  similar  to  that  of  Section  7.7.2,  a  direct  relationship  between  controller 
performance  and  parameter  quality  can  be  made.  Simply  put,  the  size  of  the  Cramer- 
Rao  bound  is  in  direct  relation  to  the  parameter  perturbation  of  the  LQR  controller, 
which  ultimately  dictates  the  necessary  control  authority  needed  to  achieve  a  desired 
level  of  control. 

The  above  discussion  reveals  the  relationship  between  the  Cramer-Rao  bound 
and  the  corresponding  magnitude  of  control  input  by  way  of  comparing  best  and  worst 
cases.  Furthermore,  it  was  shown  in  Section  7.7.2  that  the  input  frequency,  Uf,  also 
corresponds  to  the  magnitude  of  the  Cramer-Rao  bound.  Knowing  this,  it  is  then 
clear  the  Cramer-Rao  bound  and  the  magnitude  of  the  input  control  are  related  by 
way  of  the  input  frequency,  Uf.  Thus,  the  relationship  between  the  magnitude  of  the 
Cramer-Rao  bound  and  control  requirements  can  be  further  clarified  by  plotting  their 
respective  magnitudes  with  respect  to  their  input  frequencies.  This  type  of  plot  is 
presented  in  Appendix  G  for  each  rotor  blade  at  input  noise  cases  of  Sv  =  1  and  Sv  = 
1,  2,  3,  andA.  The  individual  noise  case  is  presented  to  demonstrate  the  relationship 
in  an  easy  to  read  format,  while  the  combined  noise  case  is  presented  to  demonstrate 
the  relationship  for  all  cases. 

To  demonstrate  how  the  above  described  plots  can  relate  the  magnitude  of  the 
control  input  to  the  magnitude  of  the  Cramer-Rao  bound  two  cases  will  be  considered. 
The  first  case,  as  seen  in  Figure  7.17,  will  be  for  Blade  4  where  only  the  noise  case, 
Sv  =  1,  will  be  considered.  The  second  case  will  be  for  Blade  4,  but  for  all  noise 
cases  and  is  presented  in  Figure  7.18.  The  plots  are  3D  representations  as  there  are 
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Comparison  of  Cramer  Rao  bound  and  Control  Input 
vs  Input  Frequency  for  Blade  4 


Figure  7.17:  Comparison  of  Cramer-Rao  Bound  to  Maximum  Negative  Control 

Deflection  For  Blade  4:  ojf  =  uip( 0,  0.05.  0.1,  0.2,  0.3,  0.4,  0.5),  Sv  =  1. 
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Comparison  of  Cramer  Rao  bound  and  Control  Input 
vs  Input  Frequency  for  Blade  4 


Figure  7.18:  Comparison  of  Cramer-Rao  Bound  to  Maximum  Negative  Control 

Deflection  For  Blade  4:  ojf  =  uip( 0,  0.05.0.1,  0.2,  0.3,  0.4,  0.5),  Sv  =  1,  2, 3, 4. 
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three  variables  presented.  The  plots  are  interpreted  as  a  series  of  2D  reflections  of 
the  3D  plot  with  one  2D  plot  for  each  combination  of  the  variables.  The  combination 
of  Cramer-Rao  bound  vs  input  frequency  has  already  been  discussed  in  Section  7.7.2 
and  therefore  the  corresponding  2D  plot  will  be  familiar.  Furthermore,  the  maximum 
negative  control  deflection  vs  input  frequency  was  discussed  in  Section  7.8.1  and 
the  2D  reflection  on  these  axis  will  also  be  familiar.  The  new  plot  will  be  that  of 
the  maximum  negative  control  deflection  vs  the  Cramer-Rao  bound.  It  is  this  2D 
plot  that  is  of  greatest  importance  as  it  provides  the  ability  to  directly  compare  the 
Cramer-Rao  bound  to  the  maximum  negative  control  deflection.  This,  afterall,  is 
the  ultimate  intent  of  using  the  Cramer-Rao  bound  to  evaluate  the  control  system 
performance.  By  reviewing  Case  1  in  Figure  7.17  the  Cramer-Rao  bound  vs.  Control 
deflection  reveals  that  as  the  bound  gets  large  the  corresponding  control  deflection 
gets  large.  This  is  evident  as  the  large  Cramer-Rao  bound  indicates  a  poor  parameter 
estimate,  and  with  poor  parameter  estimates  a  controller  based  upon  that  parameter 
will  therefore  have  poor  performance.  The  poor  performance,  in  this  case,  is  indicated 
by  the  controller  having  to  compensate  for  modeling  error  by  putting  in  more  control. 
This  large  control  input  is  seen  by  a  larger  negative  control  deflection.  Case  2  in 
Figure  7.18  can  be  interpreted  as  described  for  Case  1,  with  the  exception  that  Case 
2  is  for  all  noise  cases.  Case  2  reveals  that  the  control  deflection  is  generally  the  same 
despite  the  increase  in  signal  noise.  This  is  due  to  the  parameter  error,  described  by 
the  increased  Cramer-Rao  bound,  is  more  due  to  the  frequency  response  than  signal 
noise. 

By  the  results  of  the  discussions  above,  the  Cramer-Rao  bound  can  be  used  in 
control  design  by  estimating  controller  performance  based  purely  on  the  bound  of  the 
identified  parameter.  This,  in  turn,  alleviates  the  necessity  of  having  to  construct  a 
controller  at  that  test  case.  The  size  of  the  bound  threshold  can  be  determined  by 
engineering  judgement  to  determine  the  maximum  acceptable  perturbation  allowable 
in  the  system. 
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7.9  Concluding  Remarks 

This  chapter  presented  how  the  Cramer-Rao  bound  can  be  used  to  quantify  the 
performance  of  an  optimal  rotor  vibration  controller.  This  was  accomplished  by  first 
determining  the  magnitude  of  an  identified  system  parameter  variation  by  way  of  the 
direct  comparison  of  the  identified  parameter  value  to  the  corresponding  Cramer-Rao 
bound.  The  Cramer-Rao  bound  was  then  used  to  relate  how  model  perturbations 
from  the  identified  parameters  effect  controller  performance. 

In  order  to  develop  the  methods  stated  above  the  following  steps  were  first 
accomplished  in  successive  order  in  this  chapter.  For  starters,  after  briefly  stating  why 
a  vibration  controller  was  needed,  the  chapter  began  by  outlining  the  rotor  system 
LTP  equations  of  motion  based  upon  a  rigid  blade,  with  4  blades  in  total.  This  was 
followed  by  developing  the  LTP  LQR  controller  for  out  of  plane  vibration  reduction. 
The  controller  development  first  outlined  the  properties  of  a  LTP  LQR  controller 
including  the  robust  characteristics  of  guaranteed  stability  margins.  This  was  followed 
by  outlining  the  tracking  regulator  design  used  to  reduce  out  of  plane  vibrations, 
which  included  an  explanation  of  how  the  selection  of  time  periodic  gain  harmonics 
was  performed.  The  chapter  concluded  by  outlining  how  the  performance  of  each 
controller  based  upon  identified  parameters  having  increasing  Cramer-Rao  bounds 
was  quantified.  This  was  done  by  first  determining  the  magnitude  of  Cramer-Rao 
bounds  for  the  identified  system  parameter,  the  Lock  number  7.  This  was  followed 
by  relating  the  performance  of  the  controller  based  upon  the  identified  parameter  7 
to  the  corresponding  Cramer-Rao  bound.  The  next  two  paragraphs  will  detail  these 
final  two  steps  in  greater  detail,  in  the  respective  order  presented  above. 

The  validation  of  the  control  system  began  by  evaluation  the  parameters  which 
define  the  mathematical  model  in  which  it  is  synthesised  upon.  This  Section  evaluated 
how  the  variables  which  define  the  test  space,  in  this  case  input  frequency  a 7  and 
measurement  noise  Sv  affect  the  quality  of  the  system  identified  model  parameters,  in 
this  case  7.  By  reviewing  the  comparison  plots  in  Appendix  A,  it  was  shown  that  both 
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input  frequency,  Uf,  and  the  intensity  of  the  measurement  noise,  Sv,  directly  affect 
the  quality  of  the  estimate  of  the  parameter,  7.  This  was  shown  by  first  considering 
the  effect  of  the  input  frequency,  Uf,  on  the  accuracy  of  the  parameter  estimate 
by  evaluation  the  inverse  relationship  between  the  Cramer-Rao  bound  and  system 
frequency  response.  Thus,  as  frequency  response  falls  off  the  Cramer-Rao  bound 
begins  to  grow  in  magnitude.  Next,  the  effect  of  measurement  signal  noise  011  the 
magnitude  of  the  bound  was  then  evaluated.  This  was  done  to  demonstrate  the 
relationship  Sv  has  to  the  Cramer-Rao  bound  over  the  test  space  for  each  value  of 
7.  This  was  done  by  considering  how  an  increase  in  measurement  noise  reduced  the 
value  of  the  Hessian,  as  seen  in  Equation  7.57,  through  a  weighting  factor  W.  This 
ultimately  increases  the  magnitude  of  the  Cramer-Rao  bound,  as  the  square  root 
inverse  of  the  Hessian  is  the  definition  of  the  bound.  This  analysis  provided  a  method 
to  define  where  in  the  test  space  the  Cramer-Rao  bound  would  be  large,  thus  revealing 
that  identified  parameters  in  this  region  would  have  large  perturbations  from  the  true 
value. 

The  effects  of  the  parameter  accuracies  were  then  quantified  by  relating  the 
model  perturbation  they  cause  to  the  control  requirements  produced  by  the  corre¬ 
sponding  controller.  This  Section  outlined  that  while  each  controller  developed  in 
the  test  space  has  guaranteed  stability,  the  control  requirements  needed  to  offset  the 
parameter  induced  modeling  inaccuracies  may  be  excessive.  Thus,  by  relating  the 
Cramer-Rao  bound  to  the  parameter  perturbations  in  a  manner  similar  to  the  pa¬ 
rameter  validation  discussed  in  the  previous  paragraph,  a  direct  relationship  between 
controller  performance  and  parameter  quality  was  made.  This  was  done  by  noting 
that  the  size  of  the  Cramer-Rao  bound  is  in  direct  relation  to  the  parameter  perturba¬ 
tion  of  the  LQR  controller,  which  ultimately  dictates  the  necessary  control  authority 
needed  to  achieve  a  desired  level  of  control.  Thus,  the  Cramer-Rao  bound  can  be  used 
in  control  design  by  estimating  controller  performance  based  purely  on  the  bound  of 
the  identified  parameter.  This  ultimately  improves  the  process  of  control  development 
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by  alleviating  the  necessity  of  having  to  construct  a  controller  at  every  test  case  to 
determine  how  much  control  authority  the  vibration  controller  will  demand. 


This  chapter  presented  a  method  to  evaluate  the  performance  of  a  vibration  con¬ 
troller  based  upon  the  Cramer-Rao  bound.  The  next  chapter  will  present  a  summary 
of  the  overall  work  presented  in  this  document,  and  recommend  future  improvements 
to  related  works. 
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VIII.  Summary 


The  primary  objective  of  this  research  was  to  develop  a  robust  rotor  smoothing 
algorithm  based  on  system  parameters  determined  from  a  frequency  domain 
system  identification  methodology.  As  this  work  was  to  be  compatible  with  a  he¬ 
licopter  rotor  in  forward  flight,  which  is  best  represented  as  a  linear  time  periodic 
system  [12,42],  all  of  the  system  identification,  parameter  validation,  and  control 
methods  had  to  be  compatible  or  developed  for  LTP  systems.  The  performed  research 
was  extended  based  on  the  the  works  of  [12, 42]  by  introducing  a  theory  to  validate 
system  parameters  of  the  identified  model  based  on  the  Cramer- Rao  bound  [3,34], 
The  LTP  Cramer-Rao  methodology  was  then  used  to  identify  the  performance  of  the 
LTP  optimal  vibration  controller  developed  in  this  work.  This  chapter  will  review  the 
content  developed  in  each  chapter  and  make  suggestions  for  future  research. 

8.1  Summary  and  Conclusions 

This  work  began  by  outlining  the  problem  of  out  of  plane  rotor  vibrations  and 
how  the  existing  rotor  vibration  smoothing  methods  are  deficient,  as  presented  in 
chapters  1  through  111.  The  historical  methods  of  rotor  track  and  balance  and  the 
modern  rotor  smoothing  methods  were  outlined  in  terms  of  the  methods  they  employ 
to  reduce  the  out  of  plane  rotor  vibrations  caused  by  asymmetric  lift  across  the  rotor. 
The  key  aspect  of  the  review  of  these  methods  was  that  in  the  case  of  the  modern 
methods  each  failed  to  accurately  reduce  the  rotor  vibrations  due  to  an  inability  to 
accurately  model  the  rotor.  As  an  accurate  rotor  model  is  necessary  to  develop  an 
accurate  control  adjustment  solution  to  reduce  the  vibrations,  it  was  determined  that  a 
method  to  identify  a  parametric  rotor  model  from  actual  measurement,  validate  that 
model  for  accuracy,  and  then  apply  control  to  the  system  based  on  this  developed 
model  was  needed.  From  this  assessment  the  objectives  for  this  work  were  developed. 
These  objectives  were  to  model  the  rotor  system  as  a  linear  time  periodic  system  in 
state  space  form,  develop  the  Cramer-Rao  bound  for  a  LTP  system  in  state  space  form 
for  the  purpose  of  parameter  validation,  and  finally  develop  a  LTP  optimal  control 
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solution  which  uses  system  identified  parameters  to  reduce  the  out  of  plane  vibrations 
in  an  accurate  and  timely  manner. 

The  basis  of  the  work  was  the  linear  time  periodic  system.  As  such,  chapter  IV 
detailed  the  mathematical  basis  for  such  a  system.  It  was  stated  that  while  the  prin¬ 
cipals  of  linear  time  invariant  models  are  well  understood,  many  of  the  mathematical 
foundations  of  linear  time  periodic  systems  are  not.  Therefore,  this  chapter  presented 
several  of  the  required  mathematical  elements  of  linear  time  periodic  systems,  as  they 
formed  the  basis  of  the  work  covered  in  later  chapters  to  develop  the  LTP  system 
model.  The  work  of  this  chapter  began  by  describing  the  properties  of  the  Fourier 
series  as  this  was  critical  to  the  development  of  the  time  periodic  theory.  Eigenvector 
and  eigenvalues  were  briefly  described,  along  with  the  theory  of  singular  value  decom¬ 
position.  This  was  used  in  the  following  chapter  on  LTP  system  theory.  The  final 
subject  covered  was  the  Toeplitz  transformation,  which  allowed  an  infinite  dimension 
Fourier  series  to  be  recast  in  a  doubly  infinite  matrix  form.  This  was  critical  to  the 
LTP  theory  in  the  following  chapter. 

As  all  of  the  system  theory  to  be  used  to  meet  the  objectives  of  this  research 
had  to  be  cast  as  a  LTP  system  in  state  space  form,  chapter  V  outlined  the  basis 
of  this  theory.  This  chapter  described  the  formation  of  the  Linear  Time  Periodic 
system  in  terms  of  the  state  space  based  upon  the  work  of  Wereley  [42] .  The  linear 
time  periodic  state  space  representation  was  essential  to  later  chapters  in  system 
parameter  validation  and  control.  As  In  the  chapters  following  this  one,  linear  time 
periodic  analogues  of  system  parameter  validation  and  optimal  control  methods  were 
build  upon  existing  state  space  based  LTI  methods.  It  is  for  this  reason  that  a  linear 
time  periodic  state  space  operator  was  developed  in  this  chapter.  The  fundamental 
work  of  this  chapter  was  the  harmonic  balance  state  space  model,  which  allowed  for 
the  creation  of  the  LTP  state  space  operator.  This  form  was  fundamental  in  the 
creation  of  the  Cramer-Rao  bound  and  optimal  controller  of  later  chapters. 
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The  cornerstone  of  this  research  was  the  development  of  the  Cramer-Rao  lower 
bound  for  a  LTP  system,  as  presented  in  chapter  VI  This  work  was  necessary  to 
accomplish  as  the  intent  of  this  research  was  to  develop  a  more  accurate  method 
of  rotor  smoothing  by  way  of  controlling  a  verihably  accurate  rotor  system  model. 
In  the  last  chapter  the  concept  of  the  linear  time  periodic  model  was  developed, 
which  provides  the  foundation  for  an  accurate  helicopter  rotor  in  forward  flight.  The 
identified  parameters  that  populate  that  LTP  model,  however,  must  be  accurate  in 
order  to  provide  a  basis  for  an  effective  controller.  This  chapter  developed  the  Cramer- 
Rao  lower  bound  for  a  linear  periodic  system  in  order  to  validate  the  identified  system 
parameters.  This  chapter  developed  by  first  outlining  the  Cramer-Rao  inequality  and 
how  it  is  related  to  the  maximum  likelihood  estimator.  The  Cramer-Rao  lower  bound 
was  then  developed  and  adapted  for  an  LTP  system  based  on  an  LTI  analog.  A 
derivation  of  the  theory  and  methodology  required  to  generate  the  Cramer-Rao  lower 
bound  for  a  specified  parameter  in  a  linear,  time  periodic  (LTP)  system  in  state 
space  form  has  been  presented.  This  development  made  possible  the  determination 
of  the  bounded  standard  deviation  of  a  system  parameter  which  has  been  estimated 
using  any  system  identification  technique.  The  Cramer-Rao  lower  bound  represents 
the  standard  deviation  based  on  using  an  optimal  estimator,  thus  providing  a  true 
measure  of  the  accuracy  of  the  estimate.  This  development  of  parameter  accuracy 
played  a  critical  part  in  the  final  chapter,  the  development  of  the  optimal  vibration 
controller. 

The  final  aspect  of  this  work  was  the  development  of  an  optimal  control  solution 
to  reduce  out  of  plane  rotor  vibrations,  as  detailed  in  chapter  VII.  This  chapter 
culminated  the  works  of  the  previous  chapter  to  produce  a  LTP  linear  quadratic 
regulator  based  upon  a  Cramer-Rao  validated  system  parameter,  in  this  case  the  Lock 
number  7.  By  using  the  Cramer-Rao  bounds  of  parameters  developed  from  the  LTP 
frequency  domain  system  identification  method  developed  by  Hwang  [12]  ,a  technique 
to  determine  the  performance  of  the  controller  was  developed.  This  chapter  first 
developed  the  LTP  equations  of  motion  for  a  four  bladed  rotor  loosely  based  upon 
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the  AH-64  for  flap  only.  This  model  was  then  used  to  develop  a  unique  vibration 
controller  based  upon  a  feedback  reference  tracking  scheme  for  a  LTP  LQR  system. 
The  chapter  culminated  by  relating  the  size  of  the  Cramer-Rao  bound  calculated  for 
the  system  identified  parameter  to  the  controller  performance  in  terms  of  control  usage 
demanded  by  the  controller.  This  was  done  by  showing  the  variation  in  the  identified 
values  of  7  as  quantified  by  the  Cramer-Rao  bound  could  then  be  directly  related  to 
the  perturbations  in  the  model  used  to  calculate  the  optimal  control  solution.  The 
perturbations  essentially  cause  additional  noise  in  the  control  system,  that  must  be 
handled  by  using  additional  control  power.  Thus,  the  Cramer-Rao  bound  allowed  for 
quantification  of  the  controller  performance  without  having  to  calculate  the  controller. 

8.2  Recommendations  for  Future  Research 

Two  recommendations  can  be  made  to  further  the  research  presented  in  this 
work.  The  first  improvement  addresses  the  the  measurement  of  states  to  feedback  in 
the  controller,  while  the  second  makes  a  recommendation  to  improve  the  vibration 
reduction  by  suggesting  a  new  reference  tracking  signal.  These  recommendations  are 
discussed  below. 

The  first  improvement  is  to  change  the  structure  of  the  reference  tracking  vibra¬ 
tion  controller  by  changing  the  internal  structure  from  a  a  linear  quadratic  regulator 
to  a  linear  quadratic  Gaussian  controller.  This  is  important  to  address  as  the  re¬ 
quirement  to  have  perfect  measurement  of  all  system  states  is  unrealistic.  For  this 
improvement,  the  Kalman  filter  should  be  adapted  to  a  linear  time  period  system  by 
way  of  the  harmonic  balance  state  space  operator.  Thus,  a  realistic  state  estimator 
would  reduce  the  demand  to  measure  blade  flapping  angles  states  required  by  the 
current  configuration  of  the  vibration  controller. 

The  second  improvement  is  to  modify  the  selection  of  a  reference  signal  used  by 
the  vibration  controller.  The  current  reference  signal  is  based  upon  an  assumption 
that  all  blades  produce  the  same  list  at  identical  flapping  angles.  By  reducing  the 
assumption  that  the  blades  are  identical,  the  summation  of  the  lift  from  each  blade  can 


151 


be  used  as  the  error  state  feedback  to  a  reference  lift  that  is  required  for  steady  state 
operation.  This  would  improve  the  applicability  of  the  vibration  controller  proposed 
in  this  research,  as  the  assumption  of  identical  blades  is  heavily  restricting. 
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Appendix  A.  Lock  Number  Validation  Plots  via  the  Cramer- Rao 

Bound 


This  appendix  holds  all  of  the  plots  generated  by  the  Lock  parameter  validation 
by  way  of  the  Cramer-Rao  bound  in  section  7.7.2.  Four  plots  are  presented  in 
total,  each  with  a  specific  value  of  Sv  specified  by  the  range  Sv  =  1, 2,  3, 4. 
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IDENTIFIED  VALUES  OF  LOCK  PARAMETER^,  AND  COORESPONDING  CRAMER-RAO  BOUND 
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Figure  A.l:  Comparison  of  Identified  Lock  Number  to  Cramer-Rao  Bound,  Sv  —  1 
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Figure  A. 2:  Comparison  of  Identified  Lock  Number  to  Cramer-Rao  Bound,  Sv  —  2 
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Cramer-Rao  Bound 


Appendix  B.  Cramer-Rao  Bound  Plots  for  an  Individual  Rotor  Blade 


This  appendix  contains  the  Cramer-Rao  plots  for  an  individual  blade  at  a  spec¬ 
ified  measurement  signal  noise  spectral  density,  Sv.  Four  plots  are  presented  in 
total,  each  with  a  specific  value  of  Sv  specified  by  the  range  Sv  =  1, 2,  3, 4. 


CRAMER-RAO  PLOT,  NOISE  VARIANCE  =  1 


Figure  B.l:  Cramer-Rao  Bound  of  Blade,  Sv  —  1 
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Cramer-Rao  Bound 


CRAMER-RAO  PLOT,  NOISE  VARIANCE  =  2 


Figure  B.2:  Cramer-Rao  Bound  of  Blade,  Sv  —  2 
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Cramer-Rao  Bound 


CRAMER-RAO  PLOT,  NOISE  VARIANCE  =  3 


Figure  B.3:  Cramer-Rao  Bound  of  Blade,  Sv  —  3 


160 


Cramer-Rao  Bound 


CRAMER-RAO  PLOT,  NOISE  VARIANCE  =  4 


Figure  B.4:  Cramer-Rao  Bound  of  Blade,  Sv  —  4 
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Appendix  C.  Vibration  Controller  Comparison  Plots,  Noise  Case 

Sv  =  1 


his  appendix  contains  the  plots  of  the  individual  flap  angles,  calculated  LQR 


_L  gains,  control  usage,  and  vibration  controller  performance  as  computed  at  spe¬ 
cific  cases  of  input  frequency  Uf  and  noise  spectral  density  Sv.  For  the  cases  pre¬ 
sented  in  this  appendix,  the  range  of  values  of  input  and  measurement  noise  are 
ujf  =  o;p(0,  0.05.  0.1,  0.2,  0.3,  0.4,  0.5)  and  Sv  —  1  are  used,  which  results  in  7  in¬ 
dividual  cases.  Each  case  will  present  four  plots;  one  representing  the  individual  flap 
angles  of  each  blade  after  control  is  applied  to  eliminate  the  asymmetric  lift,  one  for 
calculated  LQR  gains,  one  for  control  usage,  0con  for  each  blade,  and  one  depicting 
the  vibration  controller  performance  in  terms  of  matching  the  reference  input. 
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Case  1 

Controller  Design  Parameters 


Input  Frequency,  Uf 
Uf  =  0.5u ip 


Noise  spectral  density,  Sv 

S„  =  1 


Individual  Blade  Flapping  Angles 


Figure  C.l:  Individual  Flap  Angles,  /?,  for  case  Uf  =  0.5o ;p,  Sv  —  1 
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LQR  GAINS  LQR  GAINS  LQR  GAINS 
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0  0.5  1 
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Figure  C.2:  All  LQR  Gains  for  case  u>f  =  0.5u ;p,  Sv  —  1 
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Control  Input,  0  Control  Input,  0 

r  con  r  con 


Control  For  Input  1 
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Control  For  Input  2 


Control  For  Input  3 


Control  For  Input  4 


Figure  C.3:  Control  Usage  for  case  Uf  —  0.5c Sv  —  1 
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Controller  Tracking  of  Reference  Input 


Controller  Tracking  of  Reference  Input,  Zoom  In 


Figure  C.4:  Tracking  Performance  of  Vibration  Controller  for  case  u>f  =  0.5cup,  Sv 
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Case  2 


Controller  Design  Parameters 


Input  Frequency,  Uf 


Noise  spectral  density,  Sv 


ujf  —  .4t Op 


Sv  =  1 


Individual  Blade  Flapping  Angles 


Figure  C.5:  Individual  Flap  Angles,  /?,  for  case  Uf  =  0.4o;p,  Sv  —  1 
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LQR  GAINS  LQR  GAINS  LQR  GAINS 
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TIME  IN  INTEGER  PERIODS 


Figure  C.6:  All  LQR  Gains  for  case  u>f  =  0Aujp,  Sv  —  1 


168 


Control  Input,  0  Control  Input,  0 

r  con  r  con 


Control  For  Input  1 
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Control  For  Input  2 


Control  For  Input  3 


Figure  C.7:  Control  Usage 


Control  For  Input  4 


case  u>f  =  0.4cup,  Sv  —  1 
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Controller  Tracking  of  Reference  Input 


Controller  Tracking  of  Reference  Input,  Zoom  In 


Figure  C.8:  Tracking  Performance  of  Vibration  Controller  for  case  u>f  =  0Aup,  Sv 
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Case  3 


Controller  Design  Parameters 


Input  Frequency,  Uf 


Noise  spectral  density,  Sv 


ujf  —  .  3  Up 


Sv  =  1 


Individual  Blade  Flapping  Angles 


Figure  C.9:  Individual  Flap  Angles,  /?,  for  case  Uf  =  0.3u;p,  Sv  —  1 
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LQR  GAINS  LQR  GAINS  LQR  GAINS 
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Figure  C.10:  All  LQR  Gains  for  case  Uf  =  0.3u;p,  Sv  —  1 
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Control  Input,  0  Control  Input,  0 
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Figure  C.ll:  Control  Usage 


Control  For  Input  4 


case  LUf  =  0.3ujp,  Sv  =  1 
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Controller  Tracking  of  Reference  Input 


Controller  Tracking  of  Reference  Input,  Zoom  In 


Figure  C.12:  Tracking  Performance  of  Vibration  Controller  for  case  Uf 

0.3u;p,  Sv  —  1 
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Case  4 


Controller  Design  Parameters 


Input  Frequency,  Uf 


Noise  spectral  density,  Sv 


ujf  —  .2ujp 


Sv  =  1 


Individual  Blade  Flapping  Angles 


Figure  C.13:  Individual  Flap  Angles,  (3,  for  case  t Of  =  0.2 up,  Sv  =  1 
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LQR  GAINS  LQR  GAINS  LQR  GAINS 


State  9 


0  0.5  1 
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Figure  C.14:  All  LQR  Gains  for  case  Uf  =  0.2cop ,  Sv  —  1 
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Control  Input,  0  Control  Input,  0 

r  con  r  con 


Control  For  Input  1 
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Control  For  Input  2 


Control  For  Input  3  Control  For  Input  4 


Figure  C.15:  Control  Usage  for  case  a ;/  =  0.2a ;p,  Sv  =  1 
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Controller  Tracking  of  Reference  Input 


Controller  Tracking  of  Reference  Input,  Zoom  In 


Figure  C.16:  Tracking  Performance  of  Vibration  Controller  for  case  Uf 

0.2u;p,  Sv  —  1 
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Case  5 


Controller  Design  Parameters 


Input  Frequency,  Uf 


Noise  spectral  density,  Sv 


UJf  —  .It Op 


Sv  =  1 


Individual  Blade  Flapping  Angles 


Figure  C.17:  Individual  Flap  Angles,  (3,  for  case  ujf  —  O.lt op,  Sv  =  1 


179 


LQR  GAINS  LQR  GAINS  LQR  GAINS 


Figure  C.18:  All  LQR  Gains  for  case  a ;/  =  0.1u;p,  Sv  =  1 
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Control  Input,  0  Control  Input,  0 
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Control  For  Input  1  Control  For  Input  2 


Control  For  Input  3  Control  For  Input  4 


Figure  C.19:  Control  Usage  for  case  a ;/  =  0.1up,  Sv  =  1 
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Controller  Tracking  of  Reference  Input 


Controller  Tracking  of  Reference  Input,  Zoom  In 


Figure  C.20:  Tracking  Performance  of  Vibration  Controller  for  case  Uf 

0.1  ujp,  Sv  =  1 
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Case  6 

Controller  Design  Parameters 


Input  Frequency,  Uf 
Uf  =  .05idp 


Noise  spectral  density,  Sv 

Sv  =  1 


Individual  Blade  Flapping  Angles 


Figure  C.21:  Individual  Flap  Angles,  /?,  for  case  c Of  =  0.05o;p,  Sv  —  1 
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Figure  C.22:  All  LQR  Gains  for  case  u>f  =  0.05u;p,  Sv  —  1 
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Control  Input,  0  Control  Input,  0 
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Control  For  Input  3 


Figure  C.23:  Control  Usage 


Control  For  Input  4 


for  case  u>f  =  0.05cup,  Sv  =  1 
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Controller  Tracking  of  Reference  Input 


Controller  Tracking  of  Reference  Input,  Zoom  In 


Figure  C.24:  Tracking  Performance  of  Vibration  Controller  for  case  Uf 

0.05wp,  Sv  =  1 
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Case  7 


Controller  Design  Parameters 


Input  Frequency,  Uf 


Noise  spectral  density,  Sv 


ujj  —  .Oc dp 


Sv  =  1 


Individual  Blade  Flapping  Angles 


Figure  C.25:  Individual  Flap  Angles,  (3,  for  case  t Of  =  O.Ocdp,  Sv  =  1 
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LQR  GAINS  LQR  GAINS  LQR  GAINS 
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Figure  C.26:  All  LQR  Gains  for  case  Uf  =  0.0cop,  Sv  —  1 
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Figure  C.27:  Control  Usage  for  case  a ;/  =  O.OcUp,  Sv  =  1 
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Controller  Tracking  of  Reference  Input 


Controller  Tracking  of  Reference  Input,  Zoom  In 


Figure  C.28:  Tracking  Performance  of  Vibration  Controller  for  case  Uf 

O.OtUp,  Sv  —  1 


190 


Appendix  D.  Vibration  Controller  Comparison  Plots,  Noise  Case 

Sv  =  2 


his  appendix  contains  the  plots  of  the  individual  flap  angles,  calculated  LQR 


_L  gains,  control  usage,  and  vibration  controller  performance  as  computed  at  spe¬ 
cific  cases  of  input  frequency  Uf  and  noise  spectral  density  Sv.  For  the  cases  pre¬ 
sented  in  this  appendix,  the  range  of  values  of  input  and  measurement  noise  are 
c jf  =  u;p(0,  0.05.  0.1,  0.2,  0.3,  0.4,  0.5)  and  Sv  —  2  are  used,  which  results  in  7  in¬ 
dividual  cases.  Each  case  will  present  four  plots;  one  representing  the  individual  flap 
angles  of  each  blade  after  control  is  applied  to  eliminate  the  asymmetric  lift,  one  for 
calculated  LQR  gains,  one  for  control  usage,  0con  for  each  blade,  and  one  depicting 
the  vibration  controller  performance  in  terms  of  matching  the  reference  input. 
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Case  1 

Controller  Design  Parameters 


Input  Frequency,  Uf 
Uf  =  0.5u ip 


Noise  spectral  density,  Sv 

S„  =  2 


Individual  Blade  Flapping  Angles 


Figure  D.l:  Individual  Flap  Angles,  (3,  for  case  oof  =  0.5o ;p,  Sv  =  2 
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Figure  D.2:  All  LQR  Gains  for  case  Uf  =  0.5up,  Sv  —  2 
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Figure  D.3:  Control  Usage  for  case  u>f  =  0.5a ;p,  Sv  —  2 
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Combined  Blade  Flapping  Angle, 


Controller  Tracking  of  Reference  Input 


Controller  Tracking  of  Reference  Input,  Zoom  In 


Figure  D.4:  Tracking  Performance  of  Vibration  Controller  for  case  u>f  =  0.5up,  Sv 
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Case  2 


Controller  Design  Parameters 


Input  Frequency,  Uf 


Noise  spectral  density,  Sv 


ujf  —  .4t Op 


Sv  =  2 


Individual  Blade  Flapping  Angles 


Figure  D.5:  Individual  Flap  Angles,  (3,  for  case  oof  =  0.4o;p,  Sv  =  2 
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LQR  GAINS  LQR  GAINS  LQR  GAINS 
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Figure  D.6:  All  LQR  Gains  for  case  Uf  =  0Au>p,  Sv  —  2 
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Figure  D.7:  Control  Usage  for  case  Uf  =  0.4u;p,  Sv  —  2 
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Combined  Blade  Flapping  Angle, 


Controller  Tracking  of  Reference  Input 


Controller  Tracking  of  Reference  Input,  Zoom  In 


Figure  D.8:  Tracking  Performance  of  Vibration  Controller  for  case  u>f  =  0  Auip,  Sv 
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Case  3 


Controller  Design  Parameters 


Input  Frequency,  Uf 


Noise  spectral  density,  Sv 


ujf  —  .  3  Up 


Sv  =  2 


Individual  Blade  Flapping  Angles 


Figure  D.9:  Individual  Flap  Angles,  (3,  for  case  a ;/  =  0.3ojp,  Sv  =  2 
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Figure  D.10:  All  LQR  Gains  for  case  Uf  =  0.3u;p,  Sv  =  2 
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Figure  D.ll:  Control  Usage 
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case  Uf  =  0.3u;p,  Sv  —  2 
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Controller  Tracking  of  Reference  Input 


Controller  Tracking  of  Reference  Input,  Zoom  In 


Figure  D.12:  Tracking  Performance  of  Vibration  Controller  for  case  Uf 

0.3u;p,  Sv  —  2 
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Case  4 


Controller  Design  Parameters 


Input  Frequency,  Uf 


Noise  spectral  density,  Sv 


ujf  —  .2ujp 


Sv  =  2 


Individual  Blade  Flapping  Angles 


Figure  D.13:  Individual  Flap  Angles,  /?,  for  case  Uf  =  0.2up,  Sv  —  2 
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Figure  D.14:  All  LQR  Gains  for  case  Uf  =  0.2 up,  Sv  =  2 
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Figure  D.15:  Control  Usage  for  case  Uf  —  0.2 up,  Sv  —  2 
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Controller  Tracking  of  Reference  Input 


Controller  Tracking  of  Reference  Input,  Zoom  In 


Figure  D.16:  Tracking  Performance  of  Vibration  Controller  for  case  Uf 

0.2 ujp,  Sv  =  2 
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Case  5 


Controller  Design  Parameters 


Input  Frequency,  Uf 


Noise  spectral  density,  Sv 


UJf  —  .It Op 


Sv  =  2 


Individual  Blade  Flapping  Angles 


Figure  D.17:  Individual  Flap  Angles,  /?,  for  case  Uf  =  0.1  up,  Sv  —  2 
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Figure  D.18:  All  LQR  Gains  for  case  Uf  =  O.ltUp,  Sv  =  2 
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Figure  D.19:  Control  Usage  for  case  Uf  —  0.1  u>p,  Sv  —  2 
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Controller  Tracking  of  Reference  Input 


Controller  Tracking  of  Reference  Input,  Zoom  In 


Figure  D.20:  Tracking  Performance  of  Vibration  Controller  for  case  Uf 

0.1  ujp,  Sv  —  2 
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Case  6 

Controller  Design  Parameters 


Input  Frequency,  Uf 
Uf  =  .05idp 


Noise  spectral  density,  Sv 
Sv  =  2 


Individual  Blade  Flapping  Angles 


Figure  D.21:  Individual  Flap  Angles,  (3,  for  case  a ;/  =  0.05a;p,  Sv  =  2 
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Figure  D.22:  All  LQR  Gains  for  case  Uf  =  0.05u;p,  Sv  =  2 
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Figure  D.23:  Control  Usage  for  case  Uf  =  0.05a;p,  Sv  —  2 
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Controller  Tracking  of  Reference  Input 


Controller  Tracking  of  Reference  Input,  Zoom  In 


Figure  D.24:  Tracking  Performance  of  Vibration  Controller  for  case  Uf 

0.05u;p,  Sv  —  2 
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Case  7 


Controller  Design  Parameters 


Input  Frequency,  Uf 


Noise  spectral  density,  Sv 


ujj  —  .Oc dp 


Sv  =  2 


Individual  Blade  Flapping  Angles 


Figure  D.25:  Individual  Flap  Angles,  /?,  for  case  Uf  =  O.Ocdp,  Sv  —  2 
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Figure  D.26:  All  LQR  Gains  for  case  uif  =  O.OtUp,  Sv  =  2 
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Figure  D.27:  Control  Usage  for  case  Uf  —  O.CLtp,  Sv  —  2 
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Controller  Tracking  of  Reference  Input 


Controller  Tracking  of  Reference  Input,  Zoom  In 


Figure  D.28:  Tracking  Performance  of  Vibration  Controller  for  case  Uf 

O.OtUp,  Sv  —  2 
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Appendix  E.  Vibration  Controller  Comparison  Plots,  Noise  Case 


his  appendix  contains  the  plots  of  the  individual  flap  angles,  calculated  LQR 


_L  gains,  control  usage,  and  vibration  controller  performance  as  computed  at  spe¬ 
cific  cases  of  input  frequency  Uf  and  noise  spectral  density  Sv.  For  the  cases  pre¬ 
sented  in  this  appendix,  the  range  of  values  of  input  and  measurement  noise  are 
c jf  =  u;p(0,  0.05.  0.1,  0.2,  0.3,  0.4,  0.5)  and  Sv  —  3  are  used,  which  results  in  7  in¬ 
dividual  cases.  Each  case  will  present  four  plots;  one  representing  the  individual  flap 
angles  of  each  blade  after  control  is  applied  to  eliminate  the  asymmetric  lift,  one  for 
calculated  LQR  gains,  one  for  control  usage,  0con  for  each  blade,  and  one  depicting 
the  vibration  controller  performance  in  terms  of  matching  the  reference  input. 
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Controller  Design  Parameters 
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Individual  Blade  Flapping  Angles 


Figure  E.l:  Individual  Flap  Angles,  (3,  for  case  Uf  =  0.5up,  Sv  —  3 
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Figure  E.2:  All  LQR  Gains  for  case  Uf  =  0.5u ;p,  Sv  —  3 
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Figure  E.3:  Control  Usage  for  case  Uf  =  0.5u ;p,  Sv  —  3 
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Combined  Blade  Flapping  Angle, 
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Controller  Tracking  of  Reference  Input,  Zoom  In 


Figure  E.4:  Tracking  Performance  of  Vibration  Controller  for  case  ujf  =  0.5u;p,  Sv 
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Case  2 


Controller  Design  Parameters 


Input  Frequency,  Uf 


Noise  spectral  density,  Sv 


ujf  —  .4t Op 


Sv  =  3 


Individual  Blade  Flapping  Angles 


Figure  E.5:  Individual  Flap  Angles,  (3,  for  case  a if  =  0.4a;p,  Sv  —  3 
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Figure  E.6:  All  LQR  Gains  for  case  Uf  =  0.4ujp,  Sv  —  3 
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Figure  E.7:  Control  Usage  for  case  Uf  =  0.4u;p,  Sv  —  3 
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Combined  Blade  Flapping  Angle, 
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Controller  Tracking  of  Reference  Input,  Zoom  In 


Figure  E.8:  Tracking  Performance  of  Vibration  Controller  for  case  ujf  =  QAup,  Sv 
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Case  3 
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Figure  E.9:  Individual  Flap  Angles,  (3,  for  case  a ;/  =  0.3a>p,  Sv  —  3 
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Figure  E.10:  All  LQR  Gains  for  case  Uf  =  0.3 up,  Sv  =  3 
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Figure  E.ll:  Control  Usage  for  case  Uf  =  0.3u;p,  Sv  —  3 
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Controller  Tracking  of  Reference  Input,  Zoom  In 


Figure  E.12:  Tracking  Performance  of  Vibration  Controller  for  case  uif 

0.3u;p,  Sv  —  3 
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Case  4 


Controller  Design  Parameters 


Input  Frequency,  Uf 


Noise  spectral  density,  Sv 


ujf  —  .2ujp 


Sv  =  3 


Individual  Blade  Flapping  Angles 


Figure  E.13:  Individual  Flap  Angles,  (3,  for  case  Uf  =  0.2lup ,  Sv  —  3 
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Figure  E.14:  All  LQR  Gains  for  case  Uf  =  0.2lup,  Sv  =  3 
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Figure  E.15:  Control  Usage  for  case  Uf  =  0.2 up,  Sv  —  3 
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Controller  Tracking  of  Reference  Input 


Controller  Tracking  of  Reference  Input,  Zoom  In 


Figure  E.16:  Tracking  Performance  of  Vibration  Controller  for  case  uif 

0.2u;p,  Sv  —  3 
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Case  5 


Controller  Design  Parameters 


Input  Frequency,  Uf 


Noise  spectral  density,  Sv 


UJf  —  .It Op 


Sv  =  3 


Individual  Blade  Flapping  Angles 


Figure  E.17:  Individual  Flap  Angles,  (3,  for  case  Uf  =  0.1  up,  Sv  —  3 
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Figure  E.18:  All  LQR  Gains  for  case  Uf  =  O.ltOp,  Sv  =  3 
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Figure  E.19:  Control  Usage  for  case  Uf  =  0.1  up,  Sv  —  3 
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Controller  Tracking  of  Reference  Input 


Controller  Tracking  of  Reference  Input,  Zoom  In 


Figure  E.20:  Tracking  Performance  of  Vibration  Controller  for  case  uif 

0.1  u!p,  Sv  —  3 
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Case  6 

Controller  Design  Parameters 


Input  Frequency,  Uf 
Uf  =  .05idp 


Noise  spectral  density,  Sv 
Sv  =  3 


Individual  Blade  Flapping  Angles 


Figure  E.21:  Individual  Flap  Angles,  (3,  for  case  a if  =  0.05a;p,  Sv  —  3 
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Figure  E.22:  All  LQR  Gains  for  case  Uf  =  0.05u;p,  Sv  =  3 
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Figure  E.23:  Control  Usage  for  case  Uf  =  0.05u;,p,  Sv  —  3 
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Controller  Tracking  of  Reference  Input 


Controller  Tracking  of  Reference  Input,  Zoom  In 


Figure  E.24:  Tracking  Performance  of  Vibration  Controller  for  case  uif 

0.05u;p,  Sv  —  3 
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Case  7 


Controller  Design  Parameters 


Input  Frequency,  Uf 


Noise  spectral  density,  Sv 


ujj  —  .Oc dp 


Sv  =  3 


Individual  Blade  Flapping  Angles 


Figure  E.25:  Individual  Flap  Angles,  (3,  for  case  Uf  =  O.Ocdp,  Sv  —  3 
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Figure  E.26:  All  LQR  Gains  for  case  Uf  =  O.OcUp,  Sv  =  3 
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Figure  E.27:  Control  Usage  for  case  Uf  =  0.0 up,  Sv  —  3 
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Controller  Tracking  of  Reference  Input 


Controller  Tracking  of  Reference  Input,  Zoom  In 


Figure  E.28:  Tracking  Performance  of  Vibration  Controller  for  case  uif 

O.OtUp,  Sv  —  3 
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Appendix  F.  Vibration  Controller  Comparison  Plots,  Noise  Case 

Sv  =  4 


his  appendix  contains  the  plots  of  the  individual  flap  angles,  calculated  LQR 


_L  gains,  control  usage,  and  vibration  controller  performance  as  computed  at  spe¬ 
cific  cases  of  input  frequency  Uf  and  noise  spectral  density  Sv.  For  the  cases  pre¬ 
sented  in  this  appendix,  the  range  of  values  of  input  and  measurement  noise  are 
c jf  =  u;p(0,  0.05.  0.1,  0.2,  0.3,  0.4,  0.5)  and  Sv  —  4  are  used,  which  results  in  7  in¬ 
dividual  cases.  Each  case  will  present  four  plots;  one  representing  the  individual  flap 
angles  of  each  blade  after  control  is  applied  to  eliminate  the  asymmetric  lift,  one  for 
calculated  LQR  gains,  one  for  control  usage,  0con  for  each  blade,  and  one  depicting 
the  vibration  controller  performance  in  terms  of  matching  the  reference  input. 
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Case  1 


Controller  Design  Parameters 


Input  Frequency,  Uf 


Noise  spectral  density,  Sv 


Uf  =  .5up 


Sv  =  4 


Individual  Blade  Flapping  Angles 


Figure  F.l:  Individual  Flap  Angles,  /?,  for  case  Uf  =  0.5o;p,  Sv  —  4 
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Figure  F.2:  All  LQR  Gains  for  case  Uf  =  0.5u ;p,  Sv  —  4 
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Figure  F.3:  Control  Usage  for  case  uj  —  0.5lup,  Sv  —  4 
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Combined  Blade  Flapping  Angle, 


Controller  Tracking  of  Reference  Input 


Controller  Tracking  of  Reference  Input,  Zoom  In 


Figure  F.4:  Tracking  Performance  of  Vibration  Controller  for  case  Uf  =  0.50^,  Sv 
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Case  2 


Controller  Design  Parameters 


Input  Frequency,  Uf 


Noise  spectral  density,  Sv 


ujf  —  .4t Op 


Sv  =  4 


Individual  Blade  Flapping  Angles 


Figure  F.5:  Individual  Flap  Angles,  /?,  for  case  Uf  =  0.4o;p,  Sv  —  4 
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Figure  F.6:  All  LQR  Gains  for  case  Uf  =  0.4u;p,  Sv  =  4 
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Figure  F.7:  Control  Usage  for  case  Uf  —  0Aujp,  Sv  —  4 
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Combined  Blade  Flapping  Angle, 


Controller  Tracking  of  Reference  Input 


Controller  Tracking  of  Reference  Input,  Zoom  In 


Figure  F.8:  Tracking  Performance  of  Vibration  Controller  for  case  Uf  =  O.Tup,  Sv 
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Case  3 


Controller  Design  Parameters 


Input  Frequency,  Uf 


Noise  spectral  density,  Sv 


ujf  —  .  3  Up 


Sv  =  4 


Individual  Blade  Flapping  Angles 


Figure  F.9:  Individual  Flap  Angles,  /?,  for  case  Uf  =  0.3o;p,  Sv  —  4 
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Figure  F.10:  All  LQR  Gains  for  case  ojf  =  0.3u>p,  Sv  —  4 
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Figure  K.  11:  Control  Usage  for  case  Uf  —  0.3u;p,  Sv  —  4 
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Controller  Tracking  of  Reference  Input 


Controller  Tracking  of  Reference  Input,  Zoom  In 


Figure  F.12:  Tracking  Performance  of  Vibration  Controller  for  case  Uf 

0.3u;p,  Sv  —  4 


261 


Case  4 


Controller  Design  Parameters 


Input  Frequency,  Uf 


Noise  spectral  density,  Sv 


ujf  —  .2ujp 


Sv  =  4 


Individual  Blade  Flapping  Angles 


Figure  F.13:  Individual  Flap  Angles,  /?,  for  case  Uf  =  0.2 up,  Sv  —  4 
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LQR  GAINS  LQR  GAINS  LQR  GAINS 
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Figure  F.14:  All  LQR  Gains  for  case  ojf  =  0.2up,  Sv  —  4 
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Figure  F.15:  Control  Usage  for  case  a ;/  =  0.2u>p,  Sv  —  4 
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Controller  Tracking  of  Reference  Input 


Controller  Tracking  of  Reference  Input,  Zoom  In 


Figure  F.16:  Tracking  Performance  of  Vibration  Controller  for  case  Uf 

0.2u;p,  Sv  —  4 
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Case  5 


Controller  Design  Parameters 


Input  Frequency,  Uf 


Noise  spectral  density,  Sv 


UJf  —  .It Op 


Sv  =  4 


Individual  Blade  Flapping  Angles 


Figure  F.17:  Individual  Flap  Angles,  /?,  for  case  Uf  =  0.1  up,  Sv  —  4 
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LQR  GAINS  LQR  GAINS  LQR  GAINS 


Figure  F.18:  All  LQR  Gains  for  case  Uf  =  0.1c op,  Sv  —  4 
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Figure  F.19:  Control  Usage  for  case  Uf  —  0.1(up,  Sv  —  4 
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Controller  Tracking  of  Reference  Input 


Controller  Tracking  of  Reference  Input,  Zoom  In 


Figure  F.20:  Tracking  Performance  of  Vibration  Controller  for  case  Uf 

O.lUp,  Sv  =  4 
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Case  6 

Controller  Design  Parameters 


Input  Frequency,  Uf 
Uf  =  .05idp 


Noise  spectral  density,  Sv 
Sv  =  4 


Individual  Blade  Flapping  Angles 


Figure  F.21:  Individual  Flap  Angles,  /?,  for  case  Uf  =  0.05o;p,  Sv  —  4 
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Figure  F.22:  All  LQR  Gains  for  case  Uf  =  0.05u;p,  Sv  —  4 
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Figure  F.23:  Control  Usage  for  case  Uf  =  0.05cup,  Sv  =  4 
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Controller  Tracking  of  Reference  Input 


Controller  Tracking  of  Reference  Input,  Zoom  In 


Figure  F.24:  Tracking  Performance  of  Vibration  Controller  for  case  Uf 

0.05wp,  Sv  =  4 
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Case  7 


Controller  Design  Parameters 


Input  Frequency,  Uf 


Noise  spectral  density,  Sv 


ujj  —  .Oc dp 


Sv  =  4 


Individual  Blade  Flapping  Angles 


Figure  F.25:  Individual  Flap  Angles,  /?,  for  case  Uf  =  0.0c Sv  —  4 
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Figure  F.26:  All  LQR  Gains  for  case  ojf  =  O.OtUp,  Sv  —  4 
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Figure  F.27:  Control  Usage  for  case  a ;/  =  O.OcUp,  Sv  —  4 
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Controller  Tracking  of  Reference  Input 


Controller  Tracking  of  Reference  Input,  Zoom  In 


Figure  F.28:  Tracking  Performance  of  Vibration  Controller  for  case  Uf 

O.OtUp,  Sv  —  4 
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Appendix  G.  Cramer-Rao  Bound  Relationship  to  Maximum  Control 

Requirements 


This  appendix  holds  all  of  the  plots  generated  by  comparing  the  Cramer-Rao 
bound  to  the  maximum  negative  control  deflection  for  the  range  of  input  fre¬ 
quencies  Uf  =  Up{ 0,  0.05.  0.1,  0.2,  0.3,  0.4,  0.5).  Eight  plots  are  presented  to  de¬ 
pict  the  maximum  negatve  control  deflection  for  each  blade.  Two  plots  are  presented 
per  blade,  one  for  the  noise  case  of  Sv  and  the  other  for  the  case  of  Sv  =  1,  2,  3,4. 
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Comparison  of  Cramer  Rao  bound  and  Control  Input 
vs  Input  Frequency  for  Blade  1 


Figure  G.l:  Comparison  of  Cramer- Rao  Bound  to  Maximum  Negative  Control 

Deflection  For  Blade  1:  ojf  =  uip( 0,  0.05.  0.1,  0.2,  0.3,  0.4,  0.5),  Sv  —  1 


279 


Comparison  of  Cramer  Rao  bound  and  Control  Input 
vs  Input  Frequency  for  Blade  1 


Figure  G.2:  Comparison  of  Cramer-Rao  Bound  to  Maximum  Negative  Control 

Deflection  For  Blade  1:  ojf  =  uip( 0,  0.05.  0.1,  0.2,  0.3,  0.4,  0.5),  Sv  =  1,2, 3, 4 
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Comparison  of  Cramer  Rao  bound  and  Control  Input 
vs  Input  Frequency  for  Blade  2 


Figure  G.3:  Comparison  of  Cramer-Rao  Bound  to  Maximum  Negative  Control 

Deflection  For  Blade  2:  a if  =  u;p(0,  0.05.  0.1,  0.2,  0.3,  0.4,  0.5),  Sv  —  1 
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Comparison  of  Cramer  Rao  bound  and  Control  Input 
vs  Input  Frequency  for  Blade  1 


Figure  G.4:  Comparison  of  Cramer-Rao  Bound  to  Maximum  Negative  Control 

Deflection  For  Blade  2:  a if  =  u;p(0,  0.05.  0.1,  0.2,  0.3,  0.4,  0.5),  Sv  =  1,2, 3, 4 
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Comparison  of  Cramer  Rao  bound  and  Control  Input 
vs  Input  Frequency  for  Blade  3 


Figure  G.5:  Comparison  of  Cramer-Rao  Bound  to  Maximum  Negative  Control 

Deflection  For  Blade  3:  ojf  =  uip( 0,  0.05.  0.1,  0.2,  0.3,  0.4,  0.5),  Sv  —  1 
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Comparison  of  Cramer  Rao  bound  and  Control  Input 
vs  Input  Frequency  for  Blade  3 


Figure  G.6:  Comparison  of  Cramer-Rao  Bound  to  Maximum  Negative  Control 

Deflection  For  Blade  3:  ojf  =  uip( 0,  0.05.  0.1,  0.2,  0.3,  0.4,  0.5),  Sv  =  1,2, 3, 4 
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Comparison  of  Cramer  Rao  bound  and  Control  Input 
vs  Input  Frequency  for  Blade  4 


Figure  G.7:  Comparison  of  Cramer- Rao  Bound  to  Maximum  Negative  Control 

Deflection  For  Blade  4:  ojf  =  uip( 0,  0.05.  0.1,  0.2,  0.3,  0.4,  0.5),  Sv  —  1 
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Comparison  of  Cramer  Rao  bound  and  Control  Input 
vs  Input  Frequency  for  Blade  4 


Figure  G.8:  Comparison  of  Cramer-Rao  Bound  to  Maximum  Negative  Control 

Deflection  For  Blade  4:  ojf  =  uip( 0,  0.05.  0.1,  0.2,  0.3,  0.4,  0.5),  Sv  =  1,2, 3, 4 
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